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MATHEMATICAL -STATIST I CAL  AND  DIGITAL  COMPUTER 
ANALYSIS  OF  TIME  SERIES  DATA 

The  purpose  of  this  study  is  to  develop  a series  of 
computer  programs  for  use  in  analyzing  time  series  data 
(waves),  such  as  EEG  readings.  Programs  were  used  to 
produce  reliable  estimates  of  correlations,  spectra, 
cross  spectra,  and  partial  coherences  of  multi-channel 
random  processes.  The  software  package  was  written  to 
be  easily  adaptable  to  different  sampling  rates,  amounts 
of  data,  and  numbers  of  channels.  Provisions  for  digital 
prefiltering  of  data,  detrending,  and  smoothing  (using  a 
number  of  lag  windows)  were  also  included.  Techniques  for 
estimation  of  spectra  by  fitting  single-  and  multi-channel 
autoregressive  schemes  to  sampled  data  were  also  applied 
and  found  to  yield  results  consistent  with  the  other 
methods.  All  programs  were  written  in  FORTRAN  and  run  on 
the  USNA/DTSS  computer  system. 


PREFACE 


2 


This  study  was  undertaken  as  part  of  a Trident 
Scholar  Research  project.  It  is  the  result  of  two 
semesters  of  study  during  the  academic  year  1975-76. 
The  many  hours  which  my  advisor,  Assoc.  Prof.  John 
S.  Kalme,  spent  contributing  help  and  guidance  are 
sincerely  appreciated.  I would  also  like  to  thank 
Assoc.  Prof.  Karel  Montor  for  supplying  EEG  data,  and 
Maj . David  A.  Wright  (CAF)  for  his  assistance  in  digi- 
tizing the  data. 


r 


! 

1 


3 


TABLE  OF  CONTENTS 

ABSTRACT  1 

PREFACE  2 

CHAPTER  1 4 

Elementary  Definitions  5 

Discrete  Implementation  8 

Smoothing  11 

CHAPTER  2 13 

BIBLIOGRAPHY  31 

APPENDIX  A.  Program  Listings 32 

APPENDIA  B.  Program  Outputs  *®9 


I.  Elementary  Definitions 

A sample  space  _fl  is  the  set  of  all  possible  outcomes  of  an 
experiment.  Each  possible  outcome  (BJ  is  called  an  elementary  event. 

A ( non-el ementary)  event  A in  11  is  any  subset  of  fl,  any  collection 
of  elementary  events.  A probability  measure  P defined  on  fl  is  a 
rule  which  to  each  event  A in  fl  assigns  a real  number  P(A)  (called 
the  probability  of  a)  such  that  the  following  conditions  sure 
satisfied i 

(1)  P(0)  - 0 

(2)  P(n)  - 0 

(3)  If  A^  are  pairwise  disjoint  events,  then  P(U  A^)“5.P(A^) 

A*  A. 

A random  variable  X is  a function  defined  on  a sample  space  fl, 
which  assigns  a real  or  complex  value  to  each  elementary  event  av. 
The  expectation  £(X)  of  the  random  variable  is  defined  as 

3 X(W)  P(UJ)  provided  / / Y ( uv)  / M P (uu)  < «=> 

-ft  fx 

Both  integrations  are  performed  in  a Lebesgue  sense. 

A random  process  is  a function  of  two  variables,  t and  tv, 
where  t is  a real  number  or  integer,  and  is  an  elementary  event 
in  the  sample  space  fl.  Thus,  X(t,t*/)  is  a random  process  if  t is 
allowed  to  vary  over  a n interval,  but  X(te  ,tv)  for  a fixed  tc  is  a 
random  variable.  Usually  the  second  argument  is  omitted  when 
expressing  a random  variables  X(t,^)  is  written  as  X(t). 

A random  process  is  said  to  be  stationary  in  the  wide  sense  if 
E(X(t)'X(t  + v))  depends  only  upon  v,  not  t.  This  expectation  is 
called  the  correlation  function  of  X and  is  denoted  by  R4(v). 

If  Y(t)  is  another  stationary  process  defined  on  the  same 
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sample  space,  and  E(X(t)*Y(t  ♦ v))  depends  only  on  v,  then  we  say 
that  X and  Y are  Jointly  stationary.  This  expectation  is  denoted  ^ 

by  RX)r(v)  and  is  called  the  cross  correlation  of  the  two  randon 
processes  X and  Y. 

The  Fourier  integral  of  the  autocorrelation  function  of  X is 
called  the  power  spectral  density  of  X,  and  is  evaluated  by  this 
expression: 

$,  U)  = L 

The  power  spectral  density  of  a process  indicates  the  amount  of 
energy  the  process  contains  in  any  frequency  Interval.  For  example, 

EEC  waveforms  have  much  of  their  energy  concentrated  near  10  Hz;  if 
such  a wave  were  passed  through  a 10  Hz  bandpass  filter  it  would 
lose  relatively  little  power.  Thus,  we  w expect  its  spectral 
density  function  to  have  a peak  about  10  Hz. 

Similarly,  the  Fourier  intrgral  of  the  cross  correlation  of 
X and  Y is  called  the  cross  spectral  density  of  X and  Y and  is 
denoted  by  S^Jf): 

^ , <c  -a  in  *'/' 

Sxr(/)  = JL  RirM  C ^ 

The  cross  spectral  density  represents  the  amount  of  power  shared 
by  the  two  processes  at  any  frequency.  For  example,  the  cross 
spectral  density  of  the  input  and  the  output  of  a bandpass  filter 
would  have  a very  high  cross  spectral  density  over  the  frequencies 
passed  by  the  filter.  » 

Of  course,  the  autocorrelation  and  cross  correlation  functions 
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can  be  recovered  from  the  spectral  density  and  cross  spectral 
density  functions,  respectively,  by  vising  the  inverse  Fourier 
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Integrals, 

If  the  processes  X and  Y axe  complex-valued,  the  definitions  of 
autocorrelation  and  cross  correlation  are  modified  to  use  the  con- 
jugate of  the  first  term: 

Rx  (v)  = E(  Xft)  ■ X(  *+*’)) 


Rxr  («■)*  E ( XU)  YU  tat)} 


The  coherence  between  X and  Y is  a normalized  function  of  the 
density  functions: 

y'  (y), 

“tr'r)  5,(f)-Sr(S) 

The  partial  coherence  of  X and  Y after  the  effects  of  a third  time 
series  Z have  been  removed  is  found  using  this  formidable-looking 
expression: 


5x«.  ( fi)  Syi.(S’)- 


ISrz(/')/'l 


r 


"1 


II.  Discrete  Implementation 

Suppose  we  have  N observed  values  of  two  processes  X and  Y. 
We  thus  have  X(t)  and  Y(t)  defined  only  for  tp0,l,2, . , ,N-1.  We 
then  define  the  correlation  functions  in  terns  of  an  average 
rather  than  an  expectation i 



X(X)X(*+v)  Rxr(*r)=z-5.Xfr)Y{x+*r) 

x*# 


For  most  purposes  it  is  necessary  to  determine  correlations 
only  for  values  of  v less  than  or  equal  to  a certain  limit  L. 
(Usually  L is  less  than  1Q%  or  20%  of  the  number  of  observed  data 
points  N.)  For  purposes  of  computation,  extend  the  sequences  X 
and  Y to  length  N*  (which  must  be  greater  than  N + L)  by  appending 


zeros  to  the  end;  call  these  new  sequences  X'  and  Y*.  Thus, 


X'(t)-X(t)  and  Y’(t)-Y(t)  for  t-0,l,...N-l 
but  X'(t)«  0 and  Y‘(t)»  0 for  t-N,N+l, . . .N’-l 

If  the  fast  Fourier  transform  (FFT)  is  to  be  used,  N’  must  be 
a power  of  two. 

Let  X*(u)  and  ?’(u)  be  the  Fourier  transforms  of  X'  and  Y*  t 


X'Mm  £ x '(X) c 


/N 


2 Y' <^)e 

«>« 


Now  form  a new  sequence  in  which  each  term  is  the  product  of 
the  complex  conjugate  of  the  corresponding  term  of  X*  and  the 
corresponding  term  of  Y' j then  take  the  inverse  Fourier  transform  of 
this  sequence. 
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,a .y/rt 


CW>^{x'^  r'(">e 

When  we  insert  the  above  expressions  for  X*(u)  and  Y*(u),  and  ma 
lpulate  the  (finite)  summations,  we  obtain  the  following! 

c <*>'*&[%  ft*)  rW'-'le 


»~i  *-i  "'-i 


*jr.2  2 l X'<*)  Y(+)  e 

~ Hz.0  +’* 


■».**■**■(*,  -X -Ar)/H‘ 


c.  m = 2 I Y'w  ru)[£  i e 

tm  a=#  *■  **' 


*'->  lwa  i*  (a 


It  is  easily  verified  that  the  trigonometric  expression  in  the 
brackets  is  equal  to  zero  unless  (s-t-v)  is  zero  or  an  integral 
multiple  (positive  or  negative)  of  H',  in  which  case  it  is  equal 
to  one.  Because  s and  t axe  restricted  to  the  range  0 to  N'-l, 
only  two  sets  of  (s,t)  pairs  meet  this  criterion. 

If  s-t-v13  -N ' , then  t^N'-v+s,  and  t must  be  greater  than  N'-v. 
We  are  restricting  v to  be  less  than  L,  and  N*  exceeds  N+L,  so 
N’-v  will  be  greater  than  N.  However,  X*(t)=0  if  t is  greater  than 
or  equal  to  N,  because  X*  is  only  an  extension  of  X beyond  N-l. 
Thus,  this  set  of  (s,t)  pairs  contributes  nothing  to  the  sum. 

If  s-t-v=0,  then  sst+v  and  t is  restricted  to  the  range  0 to 
N'-v-l.  Our  expression  then  reduces  to 

co v)=  mi  y'(*+~) 

Again,  because  of  the  w?y  in  which  the  original  sequences  were 
extended,  Y'(t+v)  is  nonzero  only  if  t+v  is  less  than  N,  only  if 


t is  less  than  or  equal  to  N-v-1 . Therefore,  N-v-1  say  be  taken  as 

the  upper  Unit  of  summation  and  we  have  the  following  relation t 

C C*r)  s Z X(*)Y(X+~)  = /V  RKr  (at) 
x*  o 

Thus,  by  performing  only  three  Fourier  transforms  we  can  obtain  all 
values  of  the  correlation  function  which  interest  us  simultaneously. 
One  more  transform  produces  the  cross  spectral  density  function. 
Observe  that  by  substituting  X and  X*  for  Y abd  Y ’ the  auto- 
correlation and  power  spectral  density  functions  of  X would  have 
been  obtained.  By  using  the  fast  Fourier  transform  to  perform  the 
above  calculations,  estimates  of  spectra  may  be  very  efficiently 
generated  on  a digital  computer.  This  method  was  realized  In  the 
subroutine  CROSS,  which  can  produce  either  cross  correlations  or 
autocorrelations . 


I 
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III.  Smoothing 

When  a correlation  function  R(v)  (in  this  section  R can  be 
either  an  autocorrelation  or  a cross  correlation)  is  Fourier  trans- 
formed. to  produce  a spectral  density  function,  the  explicit  relation 
between  R and  S is 

where  h is  the  sampling  interval,  and  fe  is  the  Hyquist  frequency, 
l/2h,  which  is  the  highest  frequency  which  cam  be  unambiguously 
determined  with  a given  sampling  rate.  The  S indicates  that  this 
is  a raw  estimate  of  the  density. 

Unfortunately,  the  above  expression  is  not  a consistent  esti- 
mate in  the  sense  of  mean  square  convergence.  Its  variance  does  not 
go  to  zero  an  the  number  of  sample  points  increases,  and  a graph 
of  raw  spectral  estimates  will  oscillate  wildly  about  the  true 
values  of  spectral  density. 

To  improve  the  spectral  estimates,  it  is  necessary  to  first 
"smooth”  or  average  the  spectral  values.  A very  simple  but  effec- 
tive method  is  to  replace  each  value  of  the  spectral  estimate 
with  a weighted  average  of  the  original  and  neighboring  values: 

?(0)“0.5  S(0)  + 0.5  S(l) 

S(k)“0.25  S(k-l)  + 0.5  S(k)  + 0.25  ^(k+l)  for  k=l,2, ...L-l 

S(L)=0.5  S(L-l)  + 0.5  S(L) 

This  smoothing  method  can  be  implemented  by  applying  a "window" 
to  the  correlation  function: 

R'(v)“R(v) *D(v/l) 


where  D is  a weighing  function  defined  as 
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X>CU-)  = £ ( 1+  cos  Tic  ) > !*l  i * 

This  is  known  as  a Tukey  window.  Another  possible  window, 

which  results  in  a different  degree  of  smoothing,  is  the 
modified  Tukey  window 

- O.S" V O . ‘K,  tojTTK  . 


i 


r 
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We  consider  random  processes  1 < j £ n , 

, that  is,  for  each  t,  Xftjis  a random  variable, 
where  all  X- ft)  are  defined  on  the  same  sample  space.  We 

. J 

ass  ume  E ( X(-  (*)  = PX;X.  (t)  does  not  depend  on  s, 

* J 

where  i can  equal  j.  Then 


If  i = j > y . (t)  is  called  autocorrelation  function  of 
, and  3x.y.  C5")  is  t^ie  power  spectral  density  of  X^ 
If  i/j , f?*.x  ft)  is  called  the  cross-correlation  function 
of  X-  and  Xj  , and  3x-x.  Cfj  is  t^ie  cross-spectral  density 
function . 

Assume  E X-(t)~0  for  all  j and  t. 

«*».  5XixJ  (j) f ^-‘27rft  1?  (t)dt 

- oo  1 j 

There  exists  a random  spectral  representation  of  the 

Xj  It): 

X ; (t)  - 5 CXWtK  dZ *,(>■) 


where 


zx. »zx  ao  * r 

J Xj  T— > o®  ‘A, 


T"  -L2ir\Jr 

‘ € -C 


-c  2rrt 


X:0)Jt 

J 


L 


I 
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The  (A)  are  processes  with  orthogonal  increments  and 

j — 


E{;ja><*Vk>  •/_  f»dZX]U)}  = 


= J ^ (a)  ^a)  5*^0) 


where  i can  equal  j. 

2Lx.{x)forms  a random  spectral  measure.  Let  us  consider  the 
physical  significance  of  S>v  ft). 

Consider  a linear  time  invariant  filter  with  input  a 
stationary  process  X(t)  and  output  YCt)  = 

J~e£irt>'  H(>)dZKa). 

- 

HfM  is  called  the  transfer  function.  Then 

Rvy  (t)  = f e1' 11rtx  I H Wl1  5a«  Cfc)  A 

and 

SwW*  I 

Also  , 

A 

Hence  if  we  can  get  estimates  5XV  ft)  for  £f)  and 

5XX  of  .Sx*(£)ywe  can  get  an  estimate  H(f)  Hft): 


Hft)  - 


5xr  ft) 
5„ft)  ' 


Also  , 


H&) 


1 Syyft) 


5xx  ft) 


is  an  estimate  of  the  square  of  the  gain  I Hft)  I ■ 

Most  often  we  do  not  have  explicit  expressions  for  Hft). 


I 
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Take  for  example  a system  such  as  a ship  which  acts  like  a 
black  box.  The  waves  xw  act  as  an  input  forcing  function. 
The  ship  processes  Xft;  in  some  way  and  responds  by 
pitching  as  an  output  Y(t) 

(Sometimes  we  can  write 

= S 

_ rvo  ( 


where 


HW-f  • 

is  called  impulse  response). 

We  subject  the  ship  model  to  waves  whose  spectral  distribu- 
tion ranges  over  the  frequencies  it  will  actually  encounter 
and  see  how  the  ship  pitches.  If  the  ship  encounters  waves 
of  its  natural  pitching  frequency,  the  ship  will  pitch 
badly.  In  this  case  the  design  must  be  adjusted  to  bring 
the  natural  pitching  frequency  to  some  frequency  at  which 
the  waves  normally  encountered  have  little  energy.  The 
captain  could  also  be  warned  of  the  sea  conditions  under 
which  he  will  have  to  alter  course  or  speed  to  avoid  danger- 
our  resonant  pitching.  As  another  example  consider  an  RC 
filter : 


£.ft) 


Tc 


£„(t) 


i 
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This  is  a high  pass  filter.  It  passes  high  frequencies, 
but  attenuates  low  frequencies.  One  might  wonder  how  we 
can  estimate  spectra  for  continuous  parameter  or  analog 
signals  by  sampling  at  discrete  time  points  and  by  using 
the  digital  computer  recover  the  spectra.  In  most 
applications  5*x^)  * ® for  | -f  | W for  some  \A/. 

For  high-quality  speech 

Sxx  ($■)  t O , loo  < 40  OOO  Hz 

For  multichannel  telephony 

Sxx  C j)  t o t 300  < J < 3?00  Hi- 

For  high-quality  music 

(f)  to  , (io-lg-)«H*. 

For  EEG 

SxxCf)  ~ O , /f/  > 5~C>  . 

The  following  Sampling  Theorum  holds: 

If  a function  of  time  contains  no  frequency  components 

higher  than  W hertz,  the  time  function  can  be  completely 
specified  by  determining  the  ordinates  at  a series  of  points 
spaced  — L seconds  apart.  Reconstitution  of  the  original 
time  function,  i.e.,  the  signal  wave  form  is  possible  if  the 
sample  pulses  are  passed  through  a suitable  low  pass  filter. 
This  is  important  in  time  division  multiplex  systems. 


1 f SXx(iJ  — ° for  /f|  > W , then 


I 


V*-  ~ v/  ***  '\  2rw(t  - aS/) 

Let  the  discrete  time  series  X;(t)  / t ' Q j 1 1,± 2 j . . . 
be  obtained  by  sampling  a continuous  parameter  time  series 
'Y'(-)at  time  intervals  of  length  h'  Xj  (t)  = 'Yl(t'h) 

Similarly  define  *j « I’or  £ = O >±  £ / . 

Let 

This  sampling  can  be  obtained  by  using  a PDP8-E  mini- 
computer, which  incorporates  an  analog-to-digital  converter. 


Then 


*v/£W  sk.Xj(f)e,'2rf*',<rff 


where 


£ = -«o  d 

Thus  in  order  that  * S>  ft  V - ) 

J J 

we  must  have  5y/.  y.  =•  O for  j f~J  > 

Otherwise  we  would  get  aliasing  and  frequencies  above 
will  be  folded  back.  (One  observes  this  in  old  westerns 
where  the  wheels  rotate  backwards  when  the  stagecoach 
starts  and  slows  down).  Thus  the  Nyquist  folding  frequency 

f «■  " W = pJy  ■ 

For  EEG  we  sample  at  10  msec  intervals,  or  to  get  a 
power  of  2,  we  sample  at  or  -fife  sec.  intervals  for  a 

period  of  10  sec.  or  less. 
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Let  X(i),XW,  ....  X Cm)  be  samples 

taken  at  time  intervals  of  length  la.  Assume  detrending 
has  been  performed.  We  shall  discuss  the  estimation  of 

A 

The  sample  cross -covariance  R/.y ■(*)  °f  lag  V 

between  X ■ (•)  and  X:  ( ) is  defined  to  be 
A N-v  J 

KiX.(¥)^L|X,(t)Xj(t+V)  ; V = 0,i,2j  . . .,  (H-l) 

= Xi(t)Xj(t+v)  , V* 

J t = -y+-i 

R*iXjM  = O , nun 
R<iX.Cx)=  (~v) 

A.  J J 

The  are  computed  by  using  the  FFT. 

* 0 

The  sample  cross  - spectral  density  function  or  cross 
periodog ram  between  X^C  ) anc*  Xj  (.' ) gi-ven  by 


N-' 

= hi 


A _ |'2TTl  <k 

“R  -J-  < f < J- 

/c.  , 2h~r  - xk 


*L*j 


o* 


SXlXM)-k  Z R¥jWe 

J * 1 O 

C's.  -oo 


_ ;air££.k 
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-l*T 


2h 


N 


J,  f 5 INX  (Nr(f-*-)M  c (X)  Js 

s»H*(rrt-0<«)  'X‘ 

5a 

N— ^ 00  • 


S’x.-x/P  • *' 


But  if  for  example  KO)  is  Gaussian  we  have  approximately 

p[l  „.*.(!)  > *]  * e~  s",‘  w x 


or 


Z IxiX;  w 


xV 


2. 

2- 


V^CI^G-))  » 5^(1) 

The  periodogram  itself  is  not  a good  estimate  of  the 
spectrum.  It  is  not  a consistent  estimate  in  the  sense  of 
mean  square  convergence.  Its  variance  does  not  go  to  zero 
as  N increases.  To  get  a good  estimate  of  Sx,x/f)  we 


must  "smooth"  the  periodogram. 

Take  for  the  moment  h * 
A — 


I 

3 r 


* . (A)  — S 


~ir 


£(5xcvW)  / W(*-*)s  (*)J«  y S .a; 

-TT 

c£  5 Wfa)  * 1 

j -TT  m 

N VAtf^S^ZM)  - a.ir  j*  W^CA-v-)  d<x.  +o(t) 

VMt( %,s<»)  * ^SM/»)(fkvlWJ«) 

Let  M < N Choose  M/^(  v o .i  or  O . 2_ 


We  use  estimates 


5 x.x.(f).fcZ  K(S)Rx.x.Me  f 

X‘*J  ivISM  «-AJ 

«(-*)•  *(«■)  i KCo)‘l,  $ov  |u|>± 

V**0W*))  = ■g-I-Sxi*i(f) 

I , ; k‘(VM*- 

- 1 

Approximate  confidence  intervals  for  the  *X  • GO  can 


be  found  by  using  the  fact  that  the  random  variable 


•5x;V.  (f) 


* X 


where 


XM 


r ^ pc 

We  choose  points  * — pj — 
to  estimate  the  spectra. 


**  v * 


Then  . A 

^ 'h'.m  K^m)^xcx/v) 


h:vX 

2.M 


I 

I 

I 


and  the  sums  in 
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The  computations  for 

-A 


Rx.>. 

» j 


*0  are  performed  by  efficient  use  of  the  FFT. 
We  used  the  lag  window 


KM  - 


Ms  -k 
* * * 5 i 


c-t  > 1 


Then  1 - O ■ 3 } ^ ~ ~r~  ~ 3-  "7/ 

The  program  MULSPECT  estimates  the  spectra  by  using  the 
above  lag  window. 

Another  method  I used  involves  the  use  of  a 
modified  periodogram  and  cross -periodograms  using  cosine 
taper.  Then  the  modified  periodograms  are  averaged  over 
neighboring  points. 

A modified  periodogram  is  of  the  form 


N 


U * -fr  2.  wNl(j) 


We  have  written  programs  which  involve  new  methods 
of  spectral  estimation,  namely  the  fitting  of  autoregressive 
schemes  to  given  time  series. 

The  program  AUTOREG  estimates  spectra  by  fitting  auto- 
regressive schemes  to  time  series.  We  solve  for  CLt/CLt.  . 

i Q.  (with  Cc0  -£  ) the  following  system  of 

equations 


...jA 
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^ A 

Z a3RX)<(t-s)  = o , t'-l,!, „ ^ 

5 »0 
A 


XQ)X(j*t) 


M-t 

E 


The  « are  computed  by  using  subroutine  CROSS. 

The  _ , ) are  computed  by  using  subroutine 

LEVNSN . 


Let  or1-  = Z * R xx  *) 

tr*  o 

The  spectral  estimate  is  given  by 


Sx*(f)  ■ 

Pick  3 s 2.^"  > 'yvl 
* 

Evaluate  3/x(j)  at 


K-O 


Jt 


$ = o$‘  x -n 


I 


3 3K 


for  J “ d,  *,  2,  - - - , 8 


Let 


^ •'VI  4 | - 3. 


-vn^-a. 


= - - ^ a 


2.i3-± 


- G 


5x/-iy«)-^^-a  -5,^ 


< = o 
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The  program  SPCTCLTK  estimates  spectra  by  averaging 
modified  per iodograms . 

The  program  MAl/TOREG  estimates  multichannel  spectra 
and  cross  - spectra  by  fitting  multidimensional  autoregressive 
schemes  to  the  multichannel  time  series. 

Le  t 


R,„K> 

i 

i 

i 

£ •**.(*) 


- o ) \ t j 

R(-vj  = RT(v) 


w-i 


where  5,Kf^)is  an  estimate 


Let  A(°)  ~ In*  ~ 


r 


'<  O 

o‘  - , 


Si«($) 

t 

i 

S™  (i) 

Of  SxJyK  Cy) 

7?y  -n  identity  matrix. 


1 

i 


Let  A(jy  for  j = /(  Z J . . t be  'rtX.'n. 

matrices  which  are  solutions  of  the  system  of  matrix  equations 


A A 


A(j)  RG-0  - O 


t o 

A 

Let  y 


Of  X <vv 


zero  matrix) 


ACj)  R(p 


Then  we  estimate  the  spectral  density  matrix  by 

A rJT'  A . tlirjCkl”1/'  /-'**.  A - iun'di -i 

S,(fl4f«Pe  “J  V.k  A(j)e  “J 


Let  3 =2-U  , 

A 

We  evaluate  SX(H  at  points 
for  K:  B 

A 

Let  A(j)  = O for  j s ^4Xy  . . . , 


Then 


za-i  a 


; 2T  .kt~'a 
.*  2oJK  \ / 


A(j)e  VJ  Z A(j)e 

j ~ o \L J 


All  calculations  are  performed  using  FFT. 


JPPPWW"1  IIIIRIllBjli*** 


s\  2 7 

The  elements  of  the  matrix  KM  are  computed  by 

subroutine  MAC,  outputted  in  multiplexed  form. 

The  matrices  A (1)  / A(z)  , . - /AM) 

are  computed  by  the  subroutine  MULLEV. 

Several  of  the  programs  involve  simulation  of  time 
series  with  specified  spectral  densities.  We  can  simulate 
PEG  or  any  time  series  with  almost  any  spectra.  Let 

be  a sequence  of  independent  observations  from 
N f 0 , 1 ) (white  noise). 

Let  Y(t)  = Zq,X(c^)  (digital  filter). 

Then  * - f \ 

5w(()  = |z«.e  I s«(f) 

= for  -i'-f'i 

£ i rr-wJJ.  k 


1 f we  take 


as  the  L-th 


2.  e 

■Tt-z-L 

partial  sum  of  the  Fourier  series  of  the  function 


f 


TT  ^>rr  M 

where  is  a given  function,  YC' ) will  have 

spectral  density  ^ yy(  $-  ) 

Each  X(>)  can  be  obtained  by  taking 

XM  = C?Rj  ) 

where  P, , ..  , are  random  numbers  from  the 

computer . 

We  took  Ji  - 

Simulation  was  useful  for  testing  programs. 


We  define  coherence  for  time  series  X and  Y 
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'Y-’Vn- 

fXV  * ' Sx»  (f)‘  SyyCf) 

Partial  coherence  of  XT’)  and  Y when  the  effects 
of  Z.  ('  ) are  removed 


cl  / \ 

Sxv'G)  “ szt 


| £•  /r\  I 5yzft)| 

Jrr(5)~  s«®. 


Let  v?  s be  the  effective  degrees  of  freedom, 
* x x 

Let  ~Y  (£)  be  an  estimate  of  y*  (■$•) 


t(4)  «\TPft)  , Yft)**/rvT) 


Let  to.hk~'(z)  * -ft  frl 


I-  X 


/X|  <■  i 


4 / . A 


Then  ta*h~'(  yff))  > where 

has  approximately  a normal  distribution: 

provided  "Vi  ^ ^ O H - r J “ O 

If  J * 0 , -jx  - “d  if  - ° . 

then 

/ , **<*>  F 

o-nv) ' 2'^"J 


r 


r—  2 9 

where  r -rv< , ( -m is  a random  variable  having  an  F 

distribution  with  degrees  of  freedom  "Wl,  and  -ymw 

This  can  be  used  to  test  the  null  hypothesis  HQ\ "^(^j 

against  the  alternative  ff(  : 

A 

■f  Yxrz.  ff)  has  the  same  distribu- 
tion as  '/‘xY  . 


t ion  as 


Spectra  and  cross  - spectra  and  partial  coherences  can 
be  used  to  localize  brain  tumors  and  epileptogenic  foci 
in  the  brain. 

Suppose  we  want  to  test  whether  Z drives  X and  Y. 

(Z  might  be  an  epileptogenic  focus). 

Assume 

Y *y  (*[)  , Y xz  , Y YZ  (V)  f 

Sy.ff)  , ^YY  Szz  (f)  are  significantly 

different  from  zero  over  a certain  frequency  range. 

Assume  yXZ-Y  ft)  * ° , Yyz -X  ft)  ^ but  r*Vz  4)'-°. 


Then  we  would  suspect  that  Z drives  X and  Y. 

Suppose  we  want  to  test  whether  Z drives  X, } 

Apply  the  above  analysis  to  all  possible  subsets  of  the 
recordings  taken  three  at  a time.  Suppose  all  the  spectra 

y‘z. X U)1  itji 

J 

are  nonzero  for  all  L } j ? but  Y'xx1’  Z.ft)  ~ ^ for 


Then  we  would  suspect  strongly  that  Z drives  X.  X-  V 

1 ) ■'/•V  • 

The  previously  described  partial  coherence  spectral 
analysis  can  be  extended  to  a large  number  of  data  channels 
to  test  whether  a linear  combination  of  channels  drives 
other  channels. 

The  multichannel  coherence  spectra  and  partial 
coherence  spectra  are  computed  by  the  program  SPCTBGTK. 

The  program  also  plots  the  partial  coherence  spectra  as 
well  as  the  coherence  spectra. 
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APPENDIX  A 


PROGRAM  LISTINGS 


MU  LS  PELT 

100  * DIMENSION  X(NS*LX)  ,R1  (LR*NS*NS)  ,W(M)  ,E(M),S(M*NS*NS) 

110  * DIMENSION  C(MI  *NS*NS)  ,TLAG(2*LR-U  ,Z(2*LR-I  ) ,TLAGA(LR) 

120  * RS=NUMHER  OF  CHANNELS 

130  * LX=NUMHER  OF  DATA  POINTS  FROM  EACH  CHANNEL 

140  * Ml  =M+ 1 = LENGTH  OF  TIME  LAG 

150  * LR=MAXIMUN  DESIRED  TIME  LAG  .LE.  LX 

100  * L=SMALLEST  INTEGER  SUCH  THAT  LX <=2**L 

170  * M=MAX  I MUM  LAG,  M=2**(N-1) 

180  * N IS  DEFINED  SO  THAT  2*M=2**N 
190  * H=L ENGTH  OF  SAMPLING  INTERVAL 
200  * LNXS  =LX*NS 
210  * LRNSNS  =LR*NS*NS 
220  * M1NSNS*M1 *NS*NS 

230  DIMENSION  XC  32  0, 2 ) ,R  1 ( 50, 2 , 2 ) , W ( 32  ) , F(  33 ) , S(  1 4 2 ) ,C(  33, 2 , 2 ) ,TLAG<99> 
240  *,Z< 99) ,TLAGA(50) 

250  CHARACTER  CH( 2) /" 1 " , »2»/ 

260  DATA  NS, LX, M, Ml , LR, L ,N , LXNS , LRNSNS , IN/2,320,32,33,50,9,6,640,200,2/ 
2 70  L I BR  A R Y " R EM A V" , " N LOGN " , "CROSS" , "W  PAR Z" , " M ACOR " , " COQU AD” , " MO V E” 

280  "NORMAG”  , " COHERE " , "OLDPLO"  , "GFSORT"  , " B IG"  , "SMALL"  , " PLOTl'R" 

290  H=1 ./o4. 

300  OPEN  FILE  2, "HTIOAT",  "NUMERIC" 

310  REA1)(  2 )X 

315  CALL  hPARZ(M,W) 

320  DO  1 J=1  , NS 

230  I CALL  REMAV(  LX,X(  I , J ) ) 

340  CALL  MACORINS, LX, X,LR,Rl , LXNS, LRNSNS, L) 

350  CALL  COOU  AIM  H,  NS  ,M  ,N  ,W,  R I , S,  Ml  , LR  ) 

360  CALL  COHhREUI  ,NS,S,C) 

370  DO  7 J =1  ,MI 

380  7 F(  J ) = J-1 

390  i)0  500  J=1  ,NS-1 

100  DO  _ " K=J+1,NS 

HO  HR  ITE  (0, 300)  CH(  J)  ,CH(  K) 

420  300  FORMAT  ( ' COHERENCF  FOR  CHANNELS  ',Al,'  AND  ' , A 1 ) 

130  CALL  OLOPL()(C(  I ,J,K)  ,F,M) 

150  100  FORMAT!  1 H , 2 ( F6 . 2 , 4X  , E9 . 2, 6X>  / ) 

4o0  hRITE(0,102> 

470  102  FORMAT (5 ( I H ,/)) 

180  500  CONTINUE 
400  DO  105  J=1 ,NS 
500  WRITE! 0, 1 0 7 ) CH ( J ) 

>10  107  FORMAT ! * AUTOSPECTRA  FOR  CHANNEL  /,A1  ) 

>20  CALL  OLDPLO!C( I ,J, J) ,F,M) 

540  105  WRITE ( 0, 1 02) 

550  DO  700  1 = 1 , NS-1 
560  DO  700  J=  1+1  ,NS 
■_>  70  HR  ITE  (0,00!  ) CH(  I ) ,CH(  J ) 

580  9(j  I FORMAT!  ' CROSS  CORRFLATION  BETWEEN  CHANNELS  ',A1,'  AND  ',A1) 

] >90  DO  108  L= 1 ,LR-1 
'jOO  TLAG(L)=I.-LR 

I 
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MULSPECT  (continued) 


610  108  Z(L)=RI (LR-L+I ,Jt I) 

620  LX)  100  L=LR,LR+LR-I 

6 JO  TLAG( L) =L-LR 

640  109  Z(  L) =R  I ( L-LR+1 , I ,J ) 

680  CALL  PL()TTR(Z,TLAG,LR+LR-I  ) 

6o0  WRITE (0,102) 

670  700  CONTINUE 
680  LX)  701  K=1  ,LR 
600  701  TLAGA(K) =K-I 
700  DO  801  J=1  ,NS 
710  rtRITE( 0,601 )CH( J) 

720  501  FORMAT < * AUTO  CORRELATION  FOR  CHANNEL  ',AI) 
730  CALL  PLOTTR ( R I ( I , J , J ) , TLAGA , LR ) 

740  801  HRITE(0, 102) 

750  PR  INT,  X 
760  STOP 
770  END 


3H 


5 PCTPGTK 


I oo 
1 10 
I 20 
I 40 
I 50 
100 
I 70 
IbO 
I 90 

1 05 
200 
210 

2 20 
225 
230 
240 
250 
260 
270 
280 
200 

3 00 
310 
320 
330 
340 
350 
360 
370 
380 
300 
400 

4 30 
440 
450 
460 

4 70 
4dO 
490 

5 00 
510 
520 
530 
550 
560 
570 
580 
590 
600 
620 


* 

★ 


* 


NS= EFFECT I VE 
NV»NUMBER  OF 
NB=NUMRER  OF 
NS  CANS  =N  =THE 


NUMBER  OF  SCANS  READ  IN 
CHANNELS 

FREQUENCY  RANDS  (A  POWER  OF  2) 
LEAST  POWER  OF  2 . GE.NS 


SRaSAMPLING  RATE=I/H 
X-INPUT  SER  I ES ( ARRAY) 

P-APRAY  FOR  STORING  CROSS  SPECTRA 


* ID  I M P=N  V I * ( N VI  + 1 ) * ( N R+  1 ) 

DIMENSION  X(  1024,4)  ,P(  660) 

REAL  F(38) ,42(33,4,4,4) ,SP(33) 

COMPLEX  S( 33,4,4) , SI  I 3, S 223, S I 23 , S I I 2,S332,S I 32, S22I ,S33! ,S23I 
CHARACTER  CH(4)/"I  " , "2"  , "3"  , "4"/ 


DATA  NS, NV, NR, SR, PI/800,4, 32, 64. 0,3. 14159265/ 

LIBRARY  " CCA R" 

LI  PRARY  "FAS  T1' , "TRANS",  "OLDPLO",  "GFSORT" , " BIG",  "SMALL" 


OPENF  ILE  2,  "FCHDAT"  , "NUMERIC" 

RFADC2) ( ( X< J , I ) ,J=I ,NS) , 1=1 ,NV) 

* DETREND  EACH  SERIES  RY  SUBTRACTING  FROM  EACH  SERIES 

* ITS  LEAST  SQUARES  LINEAR  REGRESSION  LINE 
FNS=NS 


THAR  =0.5*( FNS+ I . ) 

TSUMSQ=( FNS*( FNS+I .)*( FNS+FNS+I . ))/o. 
DO  76  12  = 1 , N V 


SUM=0. 


CRS PRO =0. 

DO  77  11=1 , NS 
SUM  =SUM+  X ( 1 1 , 12) 

77  CRS PRO=CRSPRO+FLOAT ( II )*X(I I ,12) 

FMEAN=SUM/FNS 

BETA  =( CRS  PRO -F  NS  *T  BA  R * F M FAN ) /( TSU  MSQ-FMS*TRAR*TBAR ) 
DO  76  11=1  , NS 

70  ■'  n , I 2)  = X(  II  ,12)  -FMEAN-BETA*(  1 1 -TBAR) 

* rilNDOH  EACH  SERIES  WITH  A COSINE  TAPER 

IR=NS/IO 

R=  IR 

DO  80  I 1 = 1 , IR 
f H = I1 

r I NT=FI 1-0.5 

a IND0N=0.5*(  I . 0-COS ( PI*FINT/R)) 

I3=NS+I  -I  I 

DO  80  12=1  , N V 

X ( I 1 , 12)  =W  INDOW*X(  11,12) 

80  X( 13, I2)=W  IND0W*X< 13, 12) 

L0G2NS  =0 
NSCANS= I 

54  IF  ( .NS . LE . NSCAN5  )G0  TO  55 
L0G2NS=L0G2NS+ I 

NSCAN S=NS CAN S+NS CANS 
GO  TO  54 

55  I F ( NS . EQ . NSCANS ) GO  TO  74 


j 


SPCTFGTK  (continued) 

630  * ZERO  OUT  THE  LAST  PART  OF  EACH  SERIES  IF  THE  NUMBER  OF  SCANS 

640  * IS  NOT  A POWER  OF  2 

650  I I REGN=NS+l 

6o0  00  7S  I 1 = II BEGN.NSCANS 

670  00  75  12=1  ,NV 

680  75  X(  II  , 12)  =0. 

690  74  CONTINUE 

700  IF ( MOD (NV,2)) 70,82,70 

710  * IF  THE  NUMBER  OF  SERIES  IS  ODD,  FILL  A DUMMY  SERIES  WITH  ZEROS 

720  70  N VI  =N  V+l 

730  00  83  I 1 =1  , N SCANS 

740  83  X C 1 1 ,N VI  ) =0.0 

750  GO  l'O  85 

760  82  N VI  =N V 

770  85  CONTINUE 

780  I D IM P=N VI  * ( N VI  + 1 ) *( NB+ I ) 

790  CALL  TRANS ( p,  ID  IM  P,  X,NSCANS,NV1  , N R,L()G2NS  ) 

800  * CROSS  SPECTRA  ESTIMATES  ARE  IN  ARRAY  P 

810  * THE  CROSS  SPECTRAL  ESTIMATES  IN  ARRAY  P ARE  SCALED  BY  MULTI PLUING 

820  * BY  Cl 

830  WNDPnR=FNS-l .25*R 

840  FSCANS=N SCANS 

850  FNR=NF 

8o0  FD=FSCANS/( FNB+FNB) 

870  C I =0. 25/ ( SR*( FD+ I . )*WND  PWR ) 

880  IR0WSP=ND+NR+2 
890  I COLS  P=(  N VI  * ( N VI  + I ) ) /2 
9 00  IS  I ZE  P=I  ROWS  P*  ICOLSP 
910  DO  95  11  = 1,  ISIZEP 
920  95  P(  I I )=CI*P(  II  ) 

925  NR  I =N P+1 
9_, . ’V)  2^0  J=l  , N V 

935  1 a=2*NBI*(NVI*(J-1  )-((  J-l  ) *(  J-2  ) ) /2-J ) 

940  DO  <;  10  K=J,NV 
945  IJK=I X+2*NB I *K 
950  DO  2D 0 I =1 , NBI 

970  200  S(  I ,J,K)=CMPLX<  P(  IJK+  1+  I- I ) , ( -I  .0)*P(  IJK+I  + I )) 

980  DO  201  J=l  , N V-l 
990  DO  201  K = J+ I ,N V 
I 000  DO  201  1=1  , NPI 

1005  CSS  = CCA B ( S ( I , J , J ) ) *CCAP<  S ( I , K , K ) ) 

1010  IF  (CSS -I  . 0E-07)  17,18,18 
1020  I 7 S( I,K,J)=(0.0,0.0) 

1030  GO  TO  201 

1040  18  S( I,K,J)=CMPLX(CCAB(S(  I , J, K > ) **2/CSS , 0.0 ) 

1050  201  CONTINUE 
1060  DO  202  Jl=l , N V-2 
I 070  DO  202  J2=J 1 + 1 , N V-l 
1080  DO  202  J3=J2+I ,NV 
1090  DO  202  1=1  ,NBI 


SPCTBGTK  (continupd) 


095  RES«REAL(S(  I,J3,J3)) 

I 00  S I I 3=S ( I , J I ,JI )*( I .0-S( I,J3,JI )) 

I 10  S223  *S( I , J2 , J2) *( 1 . -S( I , J3, J2 ) ) 

111  IF  ( RES-I  .0E-07)  901,002,902 

112  901  SI  23=S(  I,  J1  ,J2) 

113  GO  TO  903 

120  902  SI23=S(  I , J 1 ,J2)-S(I,JI  , J3  ) *C()NJG(  S ( I , J2  , J 3 ) ) /RES 
125  903  IF(CCAR(S1  I 3) *CCA P< S223 ) - I .0E-07)  601,602,602 
120  601  U2( I,JI ,J2,J3)=0.0 
127  GO  TO  202 

130  602  il  2 ( I,J1  , J 2 , J 3)  = ( CCA  R(  S 123)  **2  ) / ( CCA  P' S 113)  *CCAB(  S223  ) ) 

1 35  RES=REAL(S( I , J2,J2) ) 

140  SI  I 2=S(  I,J1  ,JI  ) *(  I . 0-S(  I,  J2 , J 1 )) 

150  S 332=S(  I,J3,J3)*(  I .-S ( I , J 3, J2 ) ) 

151  IF(  RES-I  .0E-07)  1001,1002,1002 
I 52  I 001  SI  32=S(  I,J1 ,J3) 

153  GO  TO  1003 

160  1002  SI 32=S( I,JI ,J3)-S( I , J 1 ,J2)*S( I,J2,J3)/RES 

165  1003  I F ( CCA  R ( S 1 12) *CCAB ( S 332 ) - 1 . 0E-07 ) 701,702,702 

166  701  W2(I,J1 ,J3,J2)=0.0 
1 07  GO  TO  202 

I 70  702  U2(  I , J 1 , J3,  J2  ) =(  CUAB(  SI  32  ) **2 ) /(  CCAB(  SI  I 2 ) *CCAB(  S 332 ) ) 

I 75  RES=R!-:aL(S(  I,JI  ,J1  )) 

ISO  S 22  I = S ( I , J 2 , J 2 ) * ( 1 . -S  ( I , J2 , J I ) ) 

190  S331 *S( I , J3,J3)*( 1 .-S( I ,J3,JI  )) 

191  IF( RES-I .0E-07)  1101,1102,1102 

192  1101  S23I =S( I ,J2, J3) 

193  GO  TO  I 103 

200  1102  S23I=5(I ,J2,J3)-S(  I, J I , J3> *CONJG(S ( I , J 1 , J2 ) >/RES 

201  1103  IF ( CCA  R ( S 22 1 ) * CCAB ( S 33 1 ) - I . 0E-07 ) 801,802,802 
202801  W2( I,J2,J3,JI )=0.0 

^ ' GO  TO  202 

205  802  m2  ( I,J2,J3,JI)=(  CCAB(  S23 1 ) **2  ) / ( CCAB(  S 22  1 ) *CCAR(  S33  1 )) 

210  202  CONTINUE 

213  CONI  INI  IE 

220  DO  7 J=l  ,NBI 

230  7 F ( J ) = J - 1 

240  DO  500  J = 1 ,NV-1 

250  DO  500  K=J+ 1 ,NV 

260  rtR  ITE(  0, 300 )CH(  J)  ,CH(K) 

262  300  FORMAT! 7 COHERENCE  FOR  CHANNELS  ',AI,'  AND  ',A1) 

263  DO  319  I=| , MB! 

264  SP(I)=REAL(S(I,K,J)) 

267  319  CONTINUE 

2o8  CALL  OLD PLO ( S P, F , NB) 

269  UR  ITE(0, 1 02) 

270  102  FORMAT! 25(1 H ,/)) 

280  DO  500  L=l  ,NV 

285  IF((J-L)*(K-L))41  8,500,418 
287  418  DO  512  1 = 1 ,NRl 


37 


SPCTBGTK  (continued) 


I 288  SP(  I)=^2(  I , J ,K  , L ) 

I 293  5)2  CONTINUE 

1295  WR  ITE (0,301 ) CH( J > ,CH( K) * CH(L ) 

1300  30!  FORMAT(  7 PARTIAL  COHERENCE  PET'ri  EEN  CHANNELS  ".A),'  AND 
1310  &AI,/,'  AFTER  THE  INFLUENCE  OF  CHANNEL  ',AI,'  HAS  BEEN  REMOVED'' ) 
1320  CALL  OLOPL()(SP,F,NB) 

I 330  WRITE(0,I02) 

1335  500  CONTINUE 
I 340  DO  4)7  J = 1 , N V 
I 350  W R IT  E ( 0 , 3 1 7)CH(J) 

1360  317  FORMAT ( ' AUTOSPECTRA  FOR  CHANNEL  ',A1> 

I 365  I X=2*NB I *(  N VI  *(  J-l  )-( (J-!  )*(J-2)  )/2)-l 

1370  DO  4 16  1=1 , N R I 

1390  416  SP(  I)=P(IX+I+I) 

14)0  CALL  ()LDPL()(SP(  I ) ,F,HR) 

I 420  41  7 HR  ITE(0 ,102) 

1430  STOP 
I 440  END 
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S PCTCLTK 

1 no  x NS=EFFECTI  VE  NUMBER  OF  CHANNELS  READ  IN 
IIO  * NV=NUMPER  OF  CHANNELS 

120  * N B=NUMB  ER  OF  FREQUENCY  BANDS  (A  POWER  OF  2) 

MO  * JSCANS  IS  ALSO  A POWER  OF  2 
150  * SR=SAMPL ING  RATE= I /H 
160  * X=  INPUT  SER  IES(ARRAY) 

170  * P=ARRAY  FOR  STORING  CROSS  SPECTRA 

180  DIMENSION  X(  J 02* t 2) , P< I 98) ,S( I 32) ,C( I 32) ,F( 33) 

185  DIMENSION  SI (33,2,2) 

100  CHARACTER  CH  ( 2 ) /" 1 " , "2 "/ 

2 00  EOU I VALENCE  (S,C) 

201  EQUIVALENCE  (S,SI) 

210  OAT A NS, N V, N P, JSCANS ,SR , PI/320,  2 , 32 , 1 024 , 64 ., 3 . 1 4 1 59265/ 

215  LI PRARY  "FAST" , "TRANS" , "MO VE" , "NORMAG" , "COHERE" ," PLOT" , "GFSORT" 

216  4, "PIG" , "SMALL" 

2^0  OPENFILE  2 , « NT  IDAT" , "NUMER  IC" 

230  READ!  2)  ((  X(  J,  I)  , J = 1 ,NS)  , 1=1  ,NV) 

240  * DETREND  THE  SERIES  BY  SUBTRACTING  FROM  EACH  SERIES  ITS 
250  * LEAST  SQUARES  LINEAR  REGRESSION  LINE 
2oO  FNS=NS 

270  TPAR  =0 .5* ( FNS+ 1 .0) 

280  TSUMSQ=( FNS*( FNS+ 1 .0) * ( FNS+FNS+ I .0) )/6.0 
200  DO  7o  12= I ,NV 

3 00  SUto=0 .0 
310  CRS  PR0=0 .0 
320  DO  77  11=1  , NS 
330  SUM=SU  >,+  X ( I 1 , 12) 

340  77  CRS  PRO  =CRS  PRO+FL()AT(  II  )*X(  II  ,12) 

350  F'MEAN=SUM/FNS 

360  FE'l  A =(  CRSPR0-FN5*TBAR*FMEAN)/(  TSU MSQ-FNS*TRAR*TPAR ) 

3 70  DO  78  I I =1  , NS 

_k  FREC=FMEAN+BETA*(  FLOAT ( I I ) -TBAR ) 

390  / 8 X ( 1 1 , 1 2 ) = X ( 1 1 , 1 2 ) -F R EG 

4 00  76  CONTINUE 

4 10  * WINDOW  EACH  SERIES  WITH  A COSINE  TAPER 
420  IR  =NS/ 1 0 
430  R = I R 

440  DO  79  II  =1  , IR 

450  FI! * I I 

455  F I NT = F 1 1 -0.5 

460  m I ND0H=0 .5* ( 1 . 0-C0S( PI*F INT/R ) ) 

470  I 3=NS+ I - 1 1 

480  DO  80  12=1 ,NV 

490  X<  II  , I2)=HIND0W*X(  II  ,12) 

500  80  X(  13, I2)=WIND0W*X(  13,12) 

510  79  CONTINUE 
520  L0G2NS  =0 
530  NSCANS= I 

540  54  IF (NS.LE.NSCANS)GO  TO  55 
550  LOG2NS=LOG2NS+ I 


SPCTCL'TK  (continued) 


560  NSCANS=NSCAMS+NSCANS 
570  GO  TO  54 

5B0  55  IF (NS. EQ. NSCANS ) GO  TO  74 

500  * ZERO  OUT  THE  LAST  PART  OF  EACH  SERIES  IF  THE  NUMBER  OF 

600  * SCANS  IS  NOT  A POWER  OF  2 

610  1 1 BEGN=NS+| 

o20  DO  75  1 1 = 1 1 REGN, NSCANS 

630  DO  75  12=1 ,NV 

640  75  X(  II  , I 2)  =0.0 

650  74  CONTINUE 

660  I F ( MOD( NV,2))70,82,70 

670  * IE  THE  NUMBER  OF  SERIES  IS  ODD,  FILL  A DUMMU  SERIES  WITH  ZEROS 

680  70  N V I =N  V H 

690  DO  83  11=1 , NSCANS 

700  83  X(  I I ( N VI  )=0.0 

710  GO  TO  85 

720  82  NVI =N  V 

730  85  CONTINUE 

740  I D IMP=N  VI  ★ ( h VI  + 1 ) * ( N P+ 1 ) 

750  CALL  TRANS(  P,  IDIMP.X,  NSCANS,  NVI  ,NB,L()G2NS) 

7o0  * CROSS  SPECTRUM  ESTIMATES  ARE  IN  ARRAY  P 

770  * THE  CROSS  SPECTRUM  ESTIMATES  IN  ARRAY  P ARE  SCALED  BY 

780  * MULTIPLYING  RY  C 

700  r»NDPNR=FNS-l  .25*R 

800  FSCANS=NSCAN3 

810  FN P=NB 

820  FD=FSCANS/( FNB+FNB ) 

830  Cl =0 . 25/ ( SR*( FD+ I .0)*r<NDPNR) 

840  I ROWS  P=N  R+NB+2 
850  I COLS  P=(  N VI  *(  N VI  + 1 ))/2 
860  ISIZF  P=  IRONS  P*IC()LSP 
8 -'0  95  II  =1  , ISIZEP 

880  95  P(  II  ) =C1  *P(  II  ) 

890  NRI=.IR+I 
900  DO  1 000  J=1  , NVI 
910  DO  I 000  K = J , N VI 
920  DO  I 000  1=1  , NBI 

930  SI  ( I , J,  K)  = P(  2*NBI  *(N  VI  *(  J-l  )-((J-l  ) ★(  J-2  ) ) /2+K-J)  + 1+  I-l  ) 

935  IF(J-K)99, 1000,1  000 

940  99  S I ( I,K,J)=P(2*NBI*(NVI*(J-I)-((J-I )*(J-2))/2+K-J)+I+I) 

945  1000  CONTINUE 

950  CALL  CO  H ER  E ( N P I , NVI, 5,0 

960  DO  7 J=»  , NBI 

970  7 F( J ) =J -I 

980  DO  500  J= I , N V-l 

990  DO  500  K= J+ 1 ,NV 

lOno  »tRI  TE(  0 , 300  )CH(  J ) ,CH(K) 

1 010  300  FORMAT(  ' COHERENCE  FOR  CHANNFLS  ',AI,'  AND  ',AD 
1020  CALL  PLOT ( C(  l +NBI  *(  J-I  ) +NBI  *N  VI  * ( K-l  ) ) , F , NB ) 

1040  1 00  FORMATUH  ,F5 ,0,4X,  E9.2/) 
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1050  ARITE(0,I02) 

1000  I 02  FORMAT ( 5 ( I H ,/)) 

1070  500  CONTINUE 
I 080  DO  I 05  J=1  ,NV 
1090  rfRITFIO,  I 07  ) CH < J) 

1100  107  FORMAT ( * AUTOSPECTRA  FOR  CHANNEL  A I ) 
lilO  CALL  PLOT  ( C(  I +NH I *(  J-l  )+N8 1 *N  VI  *(  J-l  ) ) , F » NB ) 
I 1 25  105  CONTINUE 
I 1 30  STOP 
I 1 40  END 
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AUl OR EG 

I 00  L IKK  ARY  "NLOGN"  , "CROSS"  , "REMA  V"  , "LE  VNSN"  , " FFT"  , "CCAR"  , "OLDPLO"  , "GFSORT" 
101  A, "PIG", "SMALL" 

110*  LR=2*NR 

120  DIMENSION  X ( 320) , R 1 ( 64 ) , A ( 64 ) , G(  33 ) , F( 33 ) 

130  COMPLEX  AC (64) 

140  * 2**LL  IS  THE  SMALLEST  POnER  OF  2 rt ITH  LX.LE.2**LL 
1 bO  DATA  LX,  NB,LR , LL, LI /320, 32 , 64 ,0 ,6/ 

160  DA'i A H/0.01  5625/ 

170  OPEN!- ILE  3,"NPUP"  ," NUMERIC" 

180  READ( 3)  ( X(  I) , XX, 1 = 1 ,320) 

190  CALL  REM A V( LX , X ) 

200  CALL  CROSS( LX , X , X, LR , R 1 , LL ) 

210  CALL  LE  VNSN ( LR , R I , A , S , M ) 

220  NR  1 =N R+ 1 

230  Ml =M+1 

2.40  DO  2 L=  I , M 1 

250  2 AC(L)=CMPLX(  A(L)  ,0.0) 

260  M2=M  1 + 1 
268  NB2=N  P+N D 
270  DO  1 1 L=M2,NP2 
280  I 1 AC ( L) =( 0. 0,0.0) 

200  * 2*NR=2**LI 
300  CALL  FFT( AC, LI  ) 

3 10  DO  3 K = 1 , N B 1 

320  3 G( K ) =( H*S ) /CABS( AC ( K ) ) **2 

330  DO  I J=1  ,33 

340  1 F( J ) =J- I 

350  CALL  OLDPLO(G,F,NR) 

370  100  FORMAT ( I HO , F5 .2, 4X , E9 . 2/ ) 

380  END 
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MAUTOREG 


I OU  LI  RRARY  " CROSS"  , "MOVE"  , "NORM AG"  , "REMA V"  , "OLD PLO" 

I 0 1 &, "GFSORT" , " R IG" , "SMALL" , "R ZERO" , "ZERO" , "COHERE" , "MA IN  V" 

I 02  A," MOVEC" , " RRA  INY" , "MATMUL" , " FFT" 

103  & , "MU LLEV" 

110*  LR=2*NR,LXNS=LX*NS 

120  * LM -SMALLEST  INTEGER  SUCH  THAT  LX.LE.2**LM 
12b  DIMENSION  PRB( 2047,2) 

1 30  DIMENSION  X(  1024,2)  ,R1  (2,2,o4)  , A ( 2, 2,64 ) , A P(  2 , 2, 64  ) , P(  2 , 2 , 64  ) 

131  &,EP( 2,2,64) , VA(2,2) , VB(2,2) , DA (2,2) ,DR(2,2) ,CA(2,2) ,CP(2,2) 

132  &,S( 33,2,2) ,C( 33,2,2) ,F( 33) 

140  COMPLEX  AC(64)  ,S  I < 2 , 2 , 33 ) ,CDP<  2, 2 ) ,ST(2,2) 

Ml  & ,C VA  ( 2 , 2 ) 

IbO  CHARACTER  CH ( 2 ) /" 1 " , "2"/ 
loO  EQUIVALENCE  (S,C) 

16b  EQUIVALENCE  <NS,NV) 

170  DATA  LX, LXNS/ 1024, 2048/ 

180  DATA  NS, NB,LR, LI ,LM/2, 32,64,6, 10/ 

190  OPENFILE  2 , " NPUB"  , "NUMER I C" 

200  READ(2)  ( PRR(  J,  I ) , J = 1 ,2047) 

210  OPENF  ILE  3,  "NCH2"  , "NUMERIC" 

220  READ(3)  ( B RR( J ,2 ) , J= I , 2047 ) 

222  DO  1 991  1=1 ,2 

224  DO  1091  J=1 ,1024 

226  1991  X(  J,  I)  = PPR(  I +2*(  J-l  ) ,1  ) 

230  DO  7 J =1 , NS 

240  7 CALL  REMA V(  LX  , X ( 1 , J ) ) 

2bO  CALL  MAC( NS, LX, X,LR,R1 ,LXNS,LM) 

2oO  CALL  MULLEVC  NS , LR, R I , A , AP, R, RP, VA , VB,DA,DP, CA, CB,M) 

2 70  DO  5 1=1  , NS 
280  DO  b J =1  ,NS 

290  b CVA ( I , J ) =CMPLX(  VA(  I , J ) ,0 .0 ) 

- ^ N PI  =N P+  I 
310  DO  1 1=1 , NS 

320  DO  1 J=1 , NS 
330  M I =M+ 1 
340  DO  2 L = 1 , Ml 

350  2 AC(l.)=CMPLX(  A(  I,  J,L)  ,0.0) 

360  M2=MI + 1 
36b  NR2=N  P+N  P 
370  DO  II  L=M2,NR2 
380  11  AC(L)=(0. 0,0.0) 

390  * 2*NB=2**LI 
400  CALL  FFT( AC, LI ) 

410  DO  I K = 1 , NB I 
420  I SI ( I , J ,K ) =AC( K ) 

430  DO  4 K=l  , NB  1 

440  CALL  MAINV(NS,SI  ( 1 ,1  , K ) , CDR) 

450  DO  3 1=1 , NS 
460  DO  3 J=l ,NS 

470  3 S 1 ( I , J , K)  =CONJG(  CD  R(  J , I ) ) 


MAUTOREG  (continued) 


480  CALL  MATMUL<NS,CDB,CVA,ST> 

490  CALL  MATMUL(  NS , ST,  SI  (I  ,1  , K)  , CDB) 

500  00  4 I = | ,NS 
510  DO  4 J=  I , NS 
520  4 SI  ( 1 ,J,K)=CDB( 1 ,J) 

5 30  DO  10  1=1  ,NS 

540  DO  10  J=  [ , NS 

550  DO  10  K=l ,NBI 

5oO  SC  K,  I , J ) =REAL( S 1 ( I , J ,K  ) ) 

570  IF(J-I)9, 10,10 

580  9 SC  K , J , I) =-A  IMAO( SI ( I, J , K ) ) 

590  10  CONTINUE 

950  CALL  COHERE ( NR  I , NS, S,C) 

960  DO  97  J=1 ,NBI 

970  97  F ( J ) = J - 1 

980  DO  500  J=1  , N V-  I 

990  DO  500  K =J+ 1 ,N  V 

1000  HR  ITE(0 , 300  )CH(  J)  ,CH(K) 

1010  300  FORMAT! * COHERENCE  FOR  CHANNELS  ',AI,'  AND  ',A1  ) 
1020  CALL  OLDPLOC  C ( I , J ,K)',  F,NB) 

1040  1 00  FORMAT! IH  , F5  .0 ,4  X,  E9 .2/) 

1050  WRITECO, 102) 

1060  102  F0RMAK25  ( I H ,/)) 

1070  5 00  CONT INUE 
1080  DO  105  J=» ,NV 
1085  WR  ITECO, 1 02) 

1090  WRITECO, 107)CH( J) 

MOO  107  FORMAT!  ' AUTOSPECTRA  FOR  CHANNEL  ',A1  ) 

1110  CALL  ()LDPL0(C(1  ,J,J),F,NB) 

1125  105  CONTINUE 
1 I 30  STO  P 
1 . • " END 

1200  SUBROUTINE  ESCALCN, A, B) 

1210  * SYMMETRIC  MATRIX  INVERSION  BY  ESCALATOR  METHOD 
1220  DIMENSION  A( N ,N ) , R( N , N) 

1 230  DO  5 1=1  ,N 
1 240  DO  5 J=l ,N 
1250  5 R( I ,J)=0.0 
1 260  B( I , I )=1  ./AC  1 ,1 ) 

1270  IFCN.EO.  I ) RETURN 
1280  DO  40  M=2,N 
1290  K=M- 1 
1300  EK=A(M,M) 

I 310  D0  10  1=  I , K 
1320  DO  10  J=l ,K 

1330  10  ED=EK-AC  M,  I ) *R(  I , J ) *A(  J,  M) 

1340  BCM,M)  =1  ./EK 
1350  DO  30  1=1  ,K 
I 360  DO  20  J=1  ,K 

1370  20  PC I , M) =B( I ,M) -d(  I,  J )*AC  J ,M ) /EK 
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I 380  30  RU,  I )=B<  I,M) 

I 390  00  40  1 = 1 ,K 
I 400  1)0  40  J=1  , K 

1410  40  B<  I, J ) =R( I , J)+B( ItM)*R(M,J)*EK 
1430  RETURN 
1430  END 

1500  SUBROUTINE  F ADDEJ(  N , A , A IN  V,DET,  ADJUG,  P) 

1510  LIRRARY  "CCAR" 

1520  DIMENSION  A ( N ,N  ) , A IN  V(  N , N ) , ADJUGC  N , N ) , P(  N ) 

1530  COMPLEX  A , A IN  V,  [)  ET,  AD  JUG,  P 
1540  COMPLEX  PN 
I 550  NN =N*N 

I 5oO  CALL  MO VEC( NN  , A , A I N V) 

15  70  DO  4 K=l  ,N 
1580  P(K)=(0. 0,0.0) 

1 590  DO  2 1=1  ,N 

1000  2 P(  K ) =P(  K ) +A  IN  V ( 1,1) 

1610  P(K)  = P(K)/FLOAT(.K) 

1620  IF ( K .EO . N ) GO  TO  5 

1630  CALL  MO VEC( N*N , A I N V, ADJUG) 

1 640  DO  3 I =1  ,N 

1650  3 ADJUG(  I,  I ) = AINV(  I,  I)-P(K) 

1 660  4 CALL  BRA  TNY(N,N, 1 , A,N,N, I , AD JUG, A INV,  1 ) 

1070  5 CALL  MOVEC(N*N, ADJUG, AINV) 

I o80  F 30  = 1 .0E-30 

1690  IF( CC A B( P( N ) ) . LT. E30) GO  TO  7 
1 700  DO  o I =1  , N 
1 710  DO  6 J=1  , N 
1720  PN  = P(  N) 

1730  o A IN  v ( I , J ) = A IN  V(  I , J ) / PN 
1 740  7 DET=P(N) 

. "o  IF  ( MOD(  N,2)  .EO.l  ) RETURN 

1 760  0FT=-DET 

1 770  DO  8 1=1  . N 

1 780  DO  8 J=l  ,N 

1790  H ADJUGC  I,J)=-AOJUC< I,J) 

1800  RETURN 
1810  END 

1900  SUBROUTINE  ESCAL1M  N , A , P,DFT) 

1910  * SYMMETRIC  MATRIX  INVERSION  BY  ESCALATOR  METHOD 
1920  * ALSO  COMPUTES  DETERMINANT  OF  A 
1930  DIMENSION  A ( N , N ) , B ( N , N ) 

1940  DO  5 I =1 ,N 
I 950  DO  5 J=l  ,N 
I960  5 R( I,J)=0.0 
1970  B(l ,1  )=l ./A ( 1 ,1  ) 

1980  DET  =A(  1,1) 

1990  I F(  N.  EO.  1 ) RETURN 
2000  DO  40  M=2,N 
2010  K=M- 1 
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2020  EK  =A  ( M , M ) 

2030  DO  10  I =1  . K 
2040  DO  10  J=l ,J 

2050  10  EK=EK-A(M,I)*R(I,J)*A(J,M) 

2060  DET=DET*EK 
2070  B ( M , M ) = I ./EK 
2080  DO  30  I =1  ,K 
2090  DO  20  J=l ,K 

2100  20  B( I ,M ) =B( I,M)-R( I,J)*A(J,M)/EK 
21  10  30  B(  M,  I ) =B(  I , M) 

21 20  DO  40  1=1 ,K 
2130  DO  40  J=l , K 

2140  40  R<  I ,J)=B(  I,J)  + R<  I , M) *B< M, J )*EK 
2150  RETURN 
2160  END 

2200  SUBROUTINE  S IMEQ( M,N, A, B, C) 

2210  * NMAX=LARG EST  VALUE  OF  N TO  BE  PROCESSED 
2220  * NONDUMMY  DIMENSION  S(NMAX.NMAX) 

2230  * FOR  EXAMPLE,  IF  NMAX=4  THEN 
2240  DIMENSION  S<4,4) 

2250  DIMENSION  A ( M , N ) ,B(N,N)  ,C(M,N) 

2260  CALL  MO VE< N*N , B ,S ) 

2270  CALL  ESCAL(N,S,P> 

2280  DO  1 1=1  ,M 
2290  DO  I J=l  ,N 
2300  A ( I , J) =0 .0 
2310  DO  1 K=l , N 

2320  I A ( I , J ) =A(  I , J)  +C(  I ,K  )*B(  K,  J) 

2330  CALL  MO VE< N*N ,S , B ) 

2340  RETURN 
2350  END 

24  f'"'  SUBROUTINE  SI MEOD(  M,N , A , R,  C,D) 

2410  * NMAX=LARGEST  VALUE  OF  N TO  RE  PROCESSED 
2420  * NCNDUMM Y DIMENSION  A(  NMAX  ,NMAX  ) 

2430  * FOR  EXAMPLE,  IF  NMAX=4,  THEN 

2440  DIMENSION  S(4,4)  * 

2450  DIMENSION  A ( M , N ) ,B(N,N),C(M,N) 

2460  CALL  MOVE(N*N, B,S) 

2470  CALL  ESCALD(N,S, R,D) 

2480  DO  I 1=1, M 
2490  DO  1 J = l , N 
2500  A(  I , J ) =0 . 0 
2510  DO  1 K=1 ,N 

2520  I A ( I , J ) =A ( I,J)+C(I,J)*B(K,J) 

2530  CALL  MO VE( N*N ,S , R) 

2540  RETURN 
2550  END 

2600  SUBROUTINE  MAC( NS ,LX, X, LR, R I , LXNS, L) 

2610  * SUBROUTINE  MAC  COMPUTES  THE  MULTICHANNEL  AUTOCORRELATION  OF 
2620  * THE  NS  CHANNEL  TIME  SERIES  X 
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I 

MAUTOREG  ( cont  i rv/^d) 

2630  * LX=LENGTH  OF  FACH  TIME  SERIES  IN  EACH  CHANNEL 
2640  * NS=NUM  RER  OF  CHANNELS 

2650  * LR«OES IRED  LENGTH  OF  CORRELATION,  LR.LE.2**L 
2660  * L IS  SMALLEST  INTEGER  SUCH  THAT  LX.LE.2**L 
26  70  * LX NS  =LX*NS 

2o80  * NON  DUMMY  DIMENSION  S(LF),  NHERE  LF.GE.ANT  ICI  PATEL)  LR 
2600  DIMENSION  X( LXNS ) ,R I ( NS, NS, LR ) 

2700  DIMENSION  S(64) 

2710  DO  I I =]  , NS 
2720  1 1 = 1 +( I-|  ) *LX 
2730  DO  I J=l , NS 
2740  Jl =l+( J-l ) *LX 

2750  CALL  CROSS ( LX , X ( I 1 ) , X ( J 1 ) , LR , S , L ) 

2760  LX)  I K =1  , LR 
2770  1 R I ( J , I , K ) =S ( K ) 

2780  RETURN 
2790  END 


ino  FUNCTION  B IG( A,M) 
110  DIMENSION  A ( M ) 

120  P=A(  I ) 

130  DO  1 K=2,M 
1 40  IF ( A( K) -B) 1,1,2 
IbO  2 P=A ( K) 

160  1 CONTINUE 
16b  B IG=P 
I 70  RETURN 
1 BO  END 


RRAINY 


100  SUBROUTINE  PRA I II  Y(  NRA  ,NCA,  LA,  A , NRR,NCB,LR,  B,  C,  LC) 

I 10  DIMENSION  A(  NRA  , NCA,  LA ) , B<  NCA , NCR,  LB)  , C(  NR  A,  NCR,  LC) 

120  COMPLEX  A,R,C 

130  * LC=LA+L  R-l 

140  CALL  l.  ERO  ( NR  A*NC R*LC,  C) 

I 50  DO  1 I =1  ,LA 
1 00  DO  I J=l ,LP 
I 70  K=  I+J-l 
1 80  DO  1 M= I , NRA 
1 90  DO  1 11  = 1 , NOB 
3 00  DO  1 L = 1 ,NCA 

210  I C(M,N,K)=C(M,N,K)+A(M,L, I)*R(L,N,J) 

220  RETURN 
230  END 


A 


L 


COHcRE 


i 

i in 
I 2 1 
I 30 
I 40 
I bO 
I 00 
I 70 
I HO 
I 90 

1 9b 
2 00 

2 0b 
210 
220 
230 
240 
250 
2oO 
270 
280 
290 
300 


SUBROUTINE  COHERE!  Ml  ,NS,S,C) 

DIMENSION  SUI  ,NS,NS)  ,C(M1  ,NS,NS) 

* EQUIVALENCE  (S,C)  IS  ALU  VI  El) 

* SUBROUTINE  COHERE  COMPUTES  THE  MAGNITUDE  AND  PHASE  ANGEL  OE 

* THE  COHERENCY,  AS  A ELL  AS  AIJTOSPECTRA,  EACH  SCALED  TO  HAVE  ITS  LARGEST 

* VALUE  UNITY 
DO  10  J p=2 , NS 
J =J  P-  I 

DO  10  K=J  P,NS 
DO  10  1=1  , Ml 

I E ( S ( I , J , J ) *S ( I , K , K ) . EO .0 . 0 ) GO  TO  10 

CO  = SORT ( A RS ( ( S ( I,J ,K)**2+S(  I , K , J ) **2 ) / ( S ( I ,J , J )*S(  I,K,K  ) ) ) ) 

I E(  A RS C S(  I , J , ,<  ) ) . LT.  I . 0E-07 ) GO  TO  101 
PH=AT AN2 ( S(  I,K,J),S(I,J,K)) 

102  C<  I,  J ,K)  =C() 

10  C( I, K,J)=180.*PH/3. 14159265 
DO  20  J = 1 , NS 

CALL  MO VE ( M 1 ,S( 1 , J, J ) ,C( 1 , J , J ) ) 

20  CALL  NORM AG( Ml , C( I , J , J ) ) 

RETURN 

101  PH  =S I GN( 1 . 5707963, S(  I,K, J ) ) 

GO  TO  102 
END 


Si 


COQUAD 


1 Oo  SUBROUTINE  COQUAD ( H , NS , M,  N , N , R I , S , M I , LR ) 

110  * SUBROUTINE  COQUAD'  COMPUTES  THE  MATRIX  OF  EMPIRICAL  AUTOSPECTRA, 

122  * COSPECTRA,  AND  QUADRATURE  SPECTRA  FROM  THE  MULTI -CHANNEL 

130  * AUTO  CORRELATION  FUNCTION 

140  * NS=NUM ^ER  OF  TIME  SERIES  OR  CHANNELS 

ISO  * to=2**( N-l  ) , M I =M+ 1 =T I ME  LENGTH  OF  CORRLAT  ION 

100  * w(  1 ) =1  , ( v,  1 ) =Q 

170  * DIMENSION  C(NM),  WHERE  NM  IS  NONDUMMY  DIMENSION  >=TW  ICE 

1 HO  * THE  MAXIMUM  LAG  M 

190  DIMENSION  ri ( M) , R 1 ( LR , NS , NS ) , S( M 1 , NS, NS ) 

2 00  COMPLEX  C(  1 00)  < 

2 10  DO  20  J=1  , NS 


220  DO  20  K =J , NS 


230  DO  10  1=1, M 

240  EVEN  =R1  ( I ,J,K)+RI ( I ,K,  J) 

250  ODD=RI ( I,J,K)-R1 ( I,K,J) 

260  R 1 ( I , K , J ) =N ( I ) *OPD 

270  10  R1  ( I,  J,K)  =.<(  I ) * EVEN 

275  20  R1 (1 , J , K) =R 1 ( 1 , J , K ) *0 .5 

2H0  DO  40  J=1  , NS 

290  DO  40  K= J , NS 

300  DO  I I =1  , M 

310  1 C(  I)=CMPLX(R!(  I,J,<)  ,RI(  I,K,JI) 

320  DO  2 I =M  I ,,M+M 
330  2 C( I) =(0. ,0.) 

340  CALL  NLOGN ( N , C , - I . , M+M ) 

350  S(  1 , J,K) =H*R EAL( C( I ) ) 

360  DO  3 1=2, Ml 

370  s S(I ,J,K)=H*(REAL(C(  I ) ) +R EALC C( M 1 +M I - I ) ) ) 

380  IF ( J -K) 7, 40, 40 
390  7 DO  4 I =2,  Ml 

4 0(  4 . T , J ) =H*(  REAL(C(  I ) ) -REALC  C(  Ml  +M  1 - I ) ) ) 

450  40  COu/INUE 

400  * S ( I , J , K ) IS  THE  COSPECTRAL  DENSITY  OF  THE  JTH  AND  KTH 
470  * CHANNEL  EVALUATED  AT  LAG  I-l  IF  K>=J,  AND  EQUAL  TO 
480  * QUADSPECTRAL  DENSITY  EVALUATED  AT  LAG  1-1  IF  K<J,  FOR 
490  RETURN 
500  END 


v 


CROSS 


I OO  SU FROUT IN E CROSS ( L HNX  , X , Y , LR  ,R  I , L) 

IIO  * 1)  IMFNS  ION  XX(.NMAX)  , YY(NMAX) 

120  X COMPLEX  CX(NMAX)  ,CY(NMAX)  ,C(  NMAX) 

130  x N.s/vx  IS  A NIONPU  ‘ >\  Y l)  I MENS  I ON  >=2  **( L+ 1 ) , MERE  2x*L.  IS  THE 
14"  * SMALLEST  POr.ER  OF  2 SUCH  THAT  LENX<=2x*L 
ISO  DIMENSION  X ( LENX)  , Y ( LENX  ) , R1  ( LR) 
loO  COMPLEX  CX(204rt) ,CY(2048) ,C(2U4d) 

16b  EQUIVALENCE  (CX,C) 

I 70  ★ LR  <=LENX 
IriO  LIBRARY  "NLOGN" 

I 8b  L2=2**( L+l  ) 

190  DO  I J = I , L EN  X 

19b  CX ( J ) =CM PLX ( X ( J ) ,0 .0) 

200  I CY(  J)=CMPLX(  Y(  J)  ,0.0) 

210  DO  2 J =LEN  X+  I ,L2 
21b  CX(J) =(0.0, 0.0) 

2^0  2 CY ( J ) =( 0.0, 0.0) 

300  CALL  NLOGN (L+l ,CX,-I ,0,L2) 

310  CAIL  NLOGN  (L+l  ,CY,-I  .0,L2) 

320  |)()  4 J-T  ,L2 

330  4 C( J)=CONJU(CX( J) )*CY( J) 

340  CALL  NLOOtK  L+l  , C, 1 ,0,L2) 

3b O DO  b J=l  , LR 

300  b R1  ( J ) =R EAL (C( J ) ) /FLOAT ( LENX) 

370  x R|  ( J)  =THE  CROSS  CORRELATION  OF  X AND  Y EVALUATED  AT 
j80  x L*G(J-I),  J = 1 , 2 , . . . M+ I 

30  0 RETURN 
4 00  END 


53 


DhlREN 


ion  SUBROUT  INF  OF  1R  Ei  <(  NS,  NV,  SCANS) 

I in*  NS=NUMRER  OF  SCANS 

ICO*  THH  FOLLOhING  IS  US  FI)  TO  DETREND  A TIME  SERIFS 
130*  PY  SUBTRACTING  ITS  LEAST  SQUARES  REGRESSION  LINE 
140*  FROM  EACH  CHANNEL  OF  AN  NV  CHANNEL  TIME  SERIES 
Ibu  DIMENSION  SCANS ( NS ,N V) 

100  r NS=NS 

170  THAR =0 ( FNS+ 1 .0) 

1 HO  TSUMSO=( ENS* (FNS+1  .0) *( 2 .0*FNS  + 1 . 0 ) ) /6 .0 
190  DO  7o  12=1 . N V 

200  SUM=0 .0 

210  CRSPR0=O.0 

220  DO  77  I l =1 , NS 

230  SUM=SUM+SCANS( I 1 , 12) 

240  CRSPRO=CHSPRO+FLOAT( I I )*SCANS< 11,12) 

2S0  77  CONTINUE 
200  FMEAN=SUM/FNS 

2 70  B ETA  =( CRS  PRO-F NS*T  RAR*FM  EAN ) / ( TSUMSQ-FNS*TBAR*TBAR ) 
280  DO  78  II  =1  , NS 

290  FREG=FMEAN+SETA*(  FLOAT ( II  ) -T RA R ) 

300  SCANS!  II  ,12)  =5 CANS ( II  , I2)-FREG 

310  78  CONTINUE 

320  76  CONTINUE 

330  RETURN 

340  END 


ESCAL 

100  Oil  R'ROll  i IN  H ESC  AL  ( N , A , B ) 

II  J * SYMMETRIC  MATRIX  INVERSION  B Y ESCALATOR  METHOD 
120  DIMENSION  A(N,N) ,R(u,N) 

I 30  00  b I =1  , N 
I 40  DO  b J=l ,N 
IbO  b «( I , J ) =0.0 
100  U(  1 , 1 ) =1  . / A ( 1,1) 

170  IE  (N. EO.  1 ) RETURN 

1 HO  00  40  A =3  , N 
190  .<=  M-  I 

200  EK  = A ( M , h\ ) 

2 10  DO  10  1=1  ,K 
22  0 00  10  J=l  ,K 

230  10  E0=  EK  - A ( M , I ) *B(  I , J ) *A ( J ,M ) 

240  ROM , M)  = I . /EK 
2bO  DO  30  1 = 1 ,K 
200  DO  20  J = 1 , K 

270  20  B ( I , M. ) =P(  I , M ) -H  ( I , J ) * A ( J , M ) / EK 
2H0  30  B ( M , I )=P( I , M ) 

200  DO  40  I =1  ,K 
300  DO  40  J=1  , K 

310  40  M( I ,J)=B(  I, J)  + P(  I ,M)*B( M, J)*EK 
320  RETURN 
330  END 


FSCALD 

inn  SUBROUTINE  HSG  AL1)(  N , A , R,  DET) 

I in  * SYMMETRIC  MATRIX  INVERSION  BY  ESCALATOR  METHOD 
120  * ALSO  COMPUTES  DETERMINANT  OF  A 
I JO  DIMENSION  A(N,N) , R(N,N> 

I 40  DO  5 1=1  , N 
ISO  DO  5 J=1 ,N 
160  b H( I , J) =0.0 
170  P(  1 , I ) = I ./ A ( I , 1 ) 

1 80  UET=A  ( 1,1) 

IRO  IF(N.F').  I ) RETURN 

2 00  DO  40  M=2,N 
210  K=M-I 

220  EK  =A ( M , M ) 

230  DO  10  1=1  , K 
240  DO  10  J= 1 , J 

250  10  EK  =EK-A  ( /<*. , I ) *B  ( I , J ) * A ( J , M ) 

260  1)ET=DE  Y*EK 
270  P(m,M)  = I./EK 
280  DO  30  1=1 ,K 
200  DO  20  J=l  , K 

300  20  B(  I,  M)  = R(  I , M)  -B(  I , J ) *A ( J , M)  / EK 
310  30  B( M, I ) =R( I , M) 

320  DO  40  I =1 , K 
330  DO  40  J=1 ,K 

340  40  R( I ,J)=B(  I,J)+R(I ,M)*B(M,J)*EK 
350  RETURN 
3oO  END 
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FADDEJ 


2*00  SUBROUTINE  FA  DDE  J (N , A , A IN V, OFT,  AD  JUG,  P) 
2410  LIBRARY  "CCAB" 

2420  I)  IMENSION  A(,.,N)  , AINV(N.N)  ,AOJUG(N,N)  , p(  N ) 
2430  COM PLi-i  X A , A IN  V , DFT , AP JUG , p 
24  40  COMPLEX  PM 
2450  NN=D*N 

24o0  CALL  MO VEC(  MM  , A ,A  IN  V) 

24  70  DO  4 K= I , N 
2480  P( K ) =( 0 . 0,0 .0) 

2490  DO  2 1=1  , II 

2500  2 P(  K)  = P(  K)+A  IN  V(  I,  I) 

2510  P(  K ) =P(  K)  /FLOAT  ( K ) 

2520  I F ( K . EO . M ) GO  TO  5 

2530  CALL  MO  VFC(  N*N , A INV , ADJUG ) 

2540  DO  3 1= 1 , N 

2550  3 APJUG(  I , I ) = A INV(  I ,1  )-P(K) 

2560  4 CALL  RR AI JY ( 3 , N , 1 , A , N , N, I , ADJUG , A I N V, 1 ) 
2570  5 CALL  MO  VEC(  N*N  , ADJUG,  A IM V) 

2580  E30= 1 . OE -JO 

25on  I F ( CCAB  ( P(  N ) ) . I.T . E30)  GO  TO  7 
2oO0  DO  0 1=1  ,N 
261 0 DO  6 J=l  ,N 
2620  PN=P(N) 

2630  6 AINV(  I,J)=A  IN \ ( I ,J)/PN 
2640  7 DET  = P(  N ) 

2650  I F ( MOD(  N , 2 ) .EO.  I ) RETURN 

2660  DFT  =-PFT 

2o70  DO  8 I=|  ,N 

2680  DO  b J= 1 , N 

2690  8 ADJUG ( I , J ) =-ADJUG ( I , J ) 

2700  RETURN 
_ / 1 0 END 


FAS  i' 


I no  SUBROUTINE  FASTI  XI  , Y I , M,N ) 

110  * SUBROUTINE  FAS']  OBTAINS  FINITE  FOUR  IERTR  ANSFORMS  OF  THE  PAIRSOF 
120  *SER I FS  X AND  Y 

130  DIMENSION  X ( I 024 ) , Y( I 024 ) , Y I ( N ) , X I ( N ) 

I 33  DO  19  J=1 ,H 

134  X(J)=XI(J) 

135  19  Y ( J ) =Y 1 ( J) 

140  N=2**M 

150  DO  1 L=  1 , M 

160  I M A X =2  ** ( 4 — L ) 

170  JPELT=2*  I MAX 

I BO  ARG=6 .283 1 853/FLOAT ( JDELT ) 

190  C =COS ( ARG) 

200  S =S I N ( ARO) 

210  U=I.O 

220  V=0.0 

230  DO  I 1=1,  IMAX 

240  DO  2 J =1 , N,  JDET.T 

250  K=J  + I MA  X 

200  XJ=X(J)+X(K) 

280  YJ=Y( J )+Y (K ) 

290  XK=X(J)-X(K) 

300  YK  =Y(  J ) -Y IK) 

310  X(  K)  =U*XK-Y*YK 

320  Y( K) =U*YK+ V*XK 

330  X(J)=XJ 

340  2 Y ( J ) =YJ 

350  UT=0*!J  -S*  V 

360  V=(J*  V+S*U 

370  1 U=UT 

380  J=1 

390  NT  =N/2 

400  I "i  A X=N- 1 

4 I 0 DO  3 I =1  , IMAX 

420  IF ( I .GE.J)  GO  TO  5 

4 30  XT=X(  J) 

440  X < J ) =X< I ) 

450  X(  I ) =XT 

460  YT=  Y ( J ) 

470  Y ( J ) =Y(  I ) 

480  Y(  I ) = YT 

490  5 K =N  1' 

5 no  4 IF(K.GE.J)  GO  TO  3 

5 10  J=J  -K 

520  =K/2 

530  GO  TO  4 

540  3 J=J+K 
542  DO  21  J=1 ,N 
544  XI (J) =X( J) 

546  21  Y 1 ( J) =Y( J) 
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FAST  (continued) 

5b  0 RETURN 
bon  end 


5? 


p 


FFT 

90  SUBROUT IMF  FFT( X,M) 

100  COMPLEX  X ( 1024),  U,  H,  T 

1 I O .><  = 2 

ICO  N V2=N/2 

I JO  N M I =N  - 1 

140  J=l 

ISO  i)()  7 1*1  f MM  I 

100  IF ( I.GE.J)  GO  TO  5 

170  T = X(J) 

I BO  X ( J ) =X  ( I ) 

I 90  X(  I ) =T 

COO  b K=MV2 

CIO  o IF(K.GE.J  ) GO  TO  7 

220  J=  J-K 

CJO  iC=  K/2 

240  GO  TO  0 

2bU  7 J=J+K 

2oO  PI=3. 1415926b 35 8979 

270  00  20  L= 1 , M 

280  LE=2**L 

200  LEI  =LE/2 

JOO  U = ( I .0,0.0) 

JIO  »'<  =CM  PLXC  COS  ( PI/FLOAT(  LEI  ))  , ( + 1 .0) *S  IN(  PI/FLO  AT  ( LEI  ) )) 

JcO  DO  20  J = l , LEI 

3J0  DO  10  I=J,N,LE 

340  I P=I  +LE  I 

350  T=X(  IP)*IJ 

360  X(IP)  = X ( I ) -T 

370  10  X(  I)=X(  D + L 

380  20  U =U*W 

300  RETURN 

400  "? 


Go 

J 


FILTER 


I nn  SU  PROU  r IN  E F I LTERt  LX  , NWT,  An , WT,  X,  LIT , Y , LN2N ) 
Jin  DIMENSION  Art(NWT) ,WT(LX) ,X(LX) ,Y(LY) 

I2U  * LY.LE.LX-NrtT+1 

13(3  * NViT=2*L+l,  2**L.i2N.  GE.LX 

I 40  130  1 J=I  , N rfi' 

ISO  I .iT(J)=AMJ) 

1 60  130  2 J=NWT+ I , L X 
I 70  2 Hi  ( J) =0.0 

I HO  CALL  CROSS ( LX , WT, X , L Y, Y.LN2N ) 

190  F N =L X 

200  DO  3 J= I , LY 

210  3 V ( J ) =F  N * Y ( J ) 

220  RETURN 
230  END 


Gl 


h. 


FTFLTtt 


I 00  SUBROUTINE  FTFLTNI  A,N«1 ,NB2, AF.M , AT.NBI ) 
I 10  DIMENSION  AC  MM  1 ) ,AT(NPI  ) 

120  CO.*';  PL  FX  A F ( N P2 ) 

130  * NP<i=2*NB=2**M.GE.Ni<1  , NB1  =NB+1 
I 40  AF(  t ) =Cf/. r*LX(  A<  1 ) ,0.0) 

130  DO  1 J=2,Nwl 

loO  1 Ar  ( J ) =CM  PL  X(  2 • *A  ( J ) ,0.0) 

I 70  DO  2 J=Nm  1 +1  ,HR2 
1M0  2 Ar( J)=(0. 0,0.0) 

1 on  CAI.L  FFT ( AF  ,3 ) 

2o o [)()  3 J=1  , N PI 

2 I 0 j AT ( J ) = R EAL ( AF ( J ) ) 

220  RETURN 

23U  END 


k ol 


UENh'LT 


i on 

1 10 

i 20 

I 30 
1 40 
I bO 
I uO 
I 70 
I dJ 

1 00 
200 

2 10 
22  0 
230 
240 
250 
200 
2 70 
280 
290 
300 
310 
320 
330 


31JWH)UTINF  <JHNi-I.T(H,A,N8l  , X,  Y,L,0,L2> 

01  MENS  ION  A ( Nd  I ) ,0(L2)  ,X(L)  , Y<  L.) 

OA i A P/0. 2831 853/ 

* t J i • I - I ri  + 1 , L 2 • 0 1: . L I 
LI =L+I 

00  b K = l ,L 

b X ( K ) =H*X(  K) 

IK)  I 1 = 2, L 

1 0(1)  =(  Y(  I)  — Y ( I -I  ))  /(  X(  I)  -X  ( I -I  )) 

0( I ) =0. 

Q( L I ) =0. 

DO  2 1=2  »NW  I 
T=2*FL0AT(  I- I ) 

TT =T *T 
A ( I ) =0 . 

DO  3 J=2 ,L I 

3 A ( I ) =A ( I ) + ( 0 ( J -1  ) -Q(  J)  )*COS(T*X(  J-l  ) )/TT 

2 A ( I ) =(  A ( I)  + ( Y(L)*SIN(T*X(L))-Y<  l)*SIN(T*X(  i)  ))/T)*2 
T=2.0*(  Y(L)*X(L) -Y(  1 )*X(  I )) 

!)()  4 J =2 , L 

4 1=3'- ( Y(  J)  -Y(  J-l  ) )*(X  ( J)  + X(  J-l  ) ) 

A ( 1 ) =T 

UPTURN 

END 


£3 


GEN  ATS 

I 00  SUBROUTINE  G ENWTS ( I D , L , i<TS , NU 2 ) 

I 10*  N*'<2=2*L+  I 

120  DIMENSION  WTS(NW2) 

130  DATA  PI/3. 14 13927/ 

Mo  IF(L.EO.O)  GO  TO  21 
I bn  r'K=L/T0-2 
160  FKK =FK+0.b 
1 70  FL=L 

I MO  DO  22  IS  = 1 ,L 
ion  S=IS 

2 no  INDEX l=(L+! )+IS 

210  HS= ( 1 .J+COSI PI *S/FL) )/(4.0*FL) 

220  22  RTS ( INDEX!  )=HS*SIN( PI*FKK*S/FL)/SIN( P I*S/ ( 2 .0*FL ) ) 

230  DO  23  IS=I  ,L 

240  INDEX  I =L+I-J3 

2b0  INDEX2=L+  I + IS 

2on  23  RTS(  INDEX  I >=RTS(  INDEX2) 

270  RTS(L+I )=FKK/FL 

2 HO  RETURN 

2°0  2 I RTS(  1 ) =1  .0 

300  RETURN 

310  END 


GFSOWT 

I nr)  SUBROUTINE  GFSORTC  G,  F , M ) 
I 1 0 DIMENSION  U(  ,F(M) 

I 23  N=F 
I 30  20  .«  = N/2 
140  IF(N) 30,40, JO 
I SO  JO  rC  = M— N 
loO  J=1 
I 70  4 I I =J 
| HO  40  L=  I+N 

i 90  I r ( G ( I ) -G ( L) )b 0,60, 60 
200  00  R=G(  I) 

210  G<  I ) =G ( L ) 

220  G(L)=B 
230  P=F(  I) 

240  r ( I ) =F ( L ) 

250  F(  L)  =H 
260  1=  I-N 

270  I h ( I - 1 )oO, 49 , 49 
260  00  J=J+1 
290  I F ( J -K) 4 1 , 4 i ,20 
300  40  RETURN 
310  END 


i 


LF  VNSN 


I 00  SUBROUTINE  LEVNSNI LR ,R , A ,Sf M) 
I 10  [)  I MENS  ION  R ( LR ) , A ( LR  ) 

I 20  V=R(  I ) 

I JO  [)=R(  2) 

140  A ( I ) = I . 0 
ISO  M=0 

100  IF  (LR. 1:0.  I ) RETURN 
170  DO  4 L=2,LR 
I MO  A ( L ) =-l)/  V 
I VO  S=V 

200  IF(  L. R0.2 ) (JO  TO  2 
210  LI =(L-2>/2 
220  L2=L I + I 
2 JO  DO  I J=2,L2 
240  HOLD  =A( J ) 

2b 0 K=L-J+I 

200  A( J)=A(J)+A(L)*A(K) 

270  I A(K)=A(K)+A( L)*HOLD 

2H0  I F ( 2*L I . EO . L-2 ) GO  TO  2 

290  A (1.2+ 1 ) * A ( L2  + I )+A(  L)  *A(  L2+1  ) 

300  2 V = V'+ A(  L)  *D 

305  ,/i=M+  I 

307  FR  INT...1 , V 

310  I F ( ( S- V)  / V -O.Ob)D,5,o 

J2  0 b IFt.V.UE.  lb) RETURN 

3J0  o IF(  L. EO .LR ) RETURN 

340  D =0 . 0 

350  DO  4 1=1  ,L 

300  t(=L-I  +2 

370  4 D=D+A( I )*R( K) 

JHO  RETURN 
39 0 b. 


I 200  SUBROUTINE  MAC(  MS , LX , X , LR , R 1 , LX  NS , L ) 

1210  * SUBROUTINE  MAC  COMPUTES  THE  MULTICHANNEL  AUTOCORRELATION  OF 
1220  * THE  NS  CHANNEL  TIME  SERIES  X 

1230  * LX=LENGI'H  OF  EACH  TIME  SERIES  IN  EACH  CHANNEL 
1240  * NS=NU,M!iER  OF  CHANNELS 

1250  x LR=DES I RED  LENGTH  OF  CORRELATION,  LR.LE.2**L 
126U  * L I,  SMALLEST  INTEGER  SUCH  THAT  LX. LE. 2**1. 

1270  * LXNS=LX*N5 

1280  * iiONDUMMY  DIMENSION  S(LF),  .-HERE  LF  .GE  . ANTI  C I RATED  LR 
I 2 DO  DIMENSION  X ( LXNS ) , R 1 ( NS , NS, LR ) 

I 3 00  DIMENSION  S ( 64 ) 

1 310  DO  I 1=1  , NS 
1320  I 1 =1  + ( I- I ) *LX 
1 330  DO  I J = l ,NS 
1 340  J I = 1 +(  J-l  )*LX 

1350  CALL  CROSS ( LX, X ( II ) ,X(JI ) ,LR,S,L) 

1 3oO  DO  I K= 1 , LR 
1 370  1 I>1  (J,  I , K)  =S(  K ) 

1380  RETURN 
1360  END 


■A  A (JOU 


IOO  SUHROU  FINE  MACOR  (NS, LX,  X, LR»  Rl  , LXNS , LRNSNS , L ) 

110  * SUPROUHNE  MACOR  COMPUTES  THE  MULTI  CHANNEL  AUTO  CORRELATION 

124  * Or  THE  NS-CHANNEL  TIME  SERIES  X 

ISO  * LX=  LENGTH  OE  EACH  LIME  SERIES  IN  EACH  CHANNEL 

14o  * NS  =NUM  PER  Or  CHANNFLS 

ISO  * LR=DES  IRED  LENGTH  OF  CORRELATION,  LR  .LE.  (LX-1) 
loo  * L IS  THE  SMALLEST  INTcGER  SUCH  THAT  LX.LE.2**L 
170  * LX NX=L X*NS , LRNSNS =LR*NS*NS 
IrsO  DIMENSION  X(  LXNS ),  R1  ( LRNSNS) 

1 00  DO  I 1=1  , NS 

200  I I =1  +(  I-l  )*LX 

21u  DO  1 J =1 , MS 

220  J I = i + ( J-l  ) *L X 

230  I J = 1 +LR* ( I -I  )+LR*NS*(  J-l  ) 

240  1 CALL  CROSS! LX, X(  I I ) , X ( Jl)  , LR  ,R  l < IJ ) ,L) 

2bO  RETURN 
2oO  END 


MA  IN  v 


I b 00  SUBROUT  I NR  1A  INI  V ( N , A , R) 

I b I 0 COM  PI . RX  A, R„DFT, ADJUG, P 
I bZO  DIMENSION  A!)JUG< 4 , 4 ) , P(  4 ) 

I b30  DIMENSION  A ( N , N ) , R(  N , N ) 

I b40  CALL  r;  AnpEJ  ( N , A , R ,DFT , ADJUG , p) 
I bbO  RETURN 
lb 60  END 


1 


G9 


MOV  EC 


100  SUBROUTINE  MOVECt LX,X, Y) 
I 10  DIMENSION  X( LX) , Y( LX ) 

MO  COMPLEX  X,Y 
I 30  DO  I 1=1  , LX 
MO  I Y(  I )=X(  I) 

150  RETURN 
loO  END 


1 


MULL tV 

2400  SUBROUTINE  MULLEV(N  ,LF,  R,  A,  AP,  B,BP,  VA  , VB,  DA.  ,D  P,  CA , CR,  M) 

2410  DIMENSION  R ( w , N ,LF) , A( N , N ,LF ) , A P( N , N ,LF ) , R(  N,N,LF) 

2420  &,BP(N,N,LF> , VA(N,N> , VB(N,N) ,DA(N,N> ,DB(N,N) ,CA(N,N) ,CB(N,N) 
2430  CALL  RZ ERO ( N*N*LF , A ) 

2440  CALL  RZERO( N*N*LF , B ) 

2450  DO  2 I =1  , N 
2460  DO  1 J=l  ,N 
2470  VA(  I ,J)=R(  I,J,I  ) 

2480  I VP(  I , J ) =R ( I ,J,1  ) 

2400  A(  I,  I,  I ) =1  . 

2500  2 P(  I,  I,  I )=l  . 

2510  CALL  ESCALDIN, VA,CR,D) 

2520  M=0 
2530  DV=D 

2540  IF(LF.EQ. 1 ) RETURN 
2550  DO  8 L=2 , LF 
2560  CALL  RZERO(N*N,DA) 

25  70  DO  5 1=1  , N 
2580  DO  4 L 1=1 ,L 
2590  LD =L-L 1+ I 
2600  DO  4 K=l  « N 
26 1 0 DO  3 J=1  , N 

2620  3 DA ( I , J ) =DA ( I , J ) —A ( I,K,LI)*R(K, J »LD) 

2630  4 CONTINUE 
2640  DO  5 J=1  ,N 
2650  5 DP(  J , I ) =DA(  I,J) 

2660  CALL  S IMEQ ( N , N , CA , VB , DA  ) 

2670  CALL  5IMEQD(N,N.CR, VA,DB,DETVA) 

2680  IF  ( (DV-DETVA)/DETVA-0.05)  1 00, 100,200 
2o90  100  IF ( M .GE  .15)  RETURN 
2700  200  D V=DETVA 
27)0  M =M  + 1 

2 720  CALL  MO  VE(  N*N*L,  A , A P) 

°730  CALL  MOVE(N*N*L,B,PP) 

2 <40  DO  7 J=l ,N 
2 750  DO  7 K=l  ,N 
2760  DO  6 L I = I , L 
2 770  LD=L-L  1+  1 
2 780  DO  6 I =1  ,N 

2700  A ( I,J,LI)=A(I,J,LI)+CA(I,K)*BP(K,J,LD) 

2800  6 R(  I,J,LI)=R(I,J,LI)+CP(  I , K)  *A  P(  K , J , LD ) 

2810  DO  7 1=1 ,N 

2820  VA( I , J ) =VA( I,J)-CA(I,K)*DB(K,J) 

2830  7 VP(I,J)=VB<  I ,J ) -CR( I ,K) *DA( K, J) 

2840  8 CONTINUE 
2850  RETURN 
2860  END 


73 


J 


NLOGN 


100  SUBROUTINE  NLOGN ( N , X,  S ION , LX  ) 

101*  NMAX=LARGEST  VALUE  OF  N TO  BE  PROCESSED 
102*  NONDUMMY  DIMENSION  M(NMAX) 

103*  FOR  EXAMPLE,  IF  NMAX=I00,  THEN 
110  DIMENSION  M(IOO) 

119*  DIMENSION  X(2**N) 

120  DIMENSION  X(LX) 

130  COMPLEX  X,  rtK,  HOLD,  Q 
140  DO  I 1=1 ,N 

1 SO  I M(  I)  =2**(N-  I) 

1 oO  DO  4 L = l , N 

170  NBL0CK=2**( L- I ) 

I 80  LBLOCK=LX/N BLOCK 

190  LBHALF=LBL0CK/2 

200  K=0 

210  DO  4 I RLOCK=l ,N BLOCK 

220  FK=K 
230  FLX=LX 

240  V=SIGN*6.2831853*FK/FLX 

250  NK=CMPLX(COS(  V)  ,SIN(  V)) 

2oO  I START=LRLOCK*( IBLOCK-I ) 

270  DO  2 1=1 , LRH ALF 

280  J=  ISTART+ I 
290  JH  =J+LBHALF 
300  0 = X(JH)*rtK 

310  X ( JH  ) =X(  J ) -0 

320  X(J)=X(J)+0 

330  2 CONTINUE 
340  DO  3 1=2, N 
350  1 1 =1 

360  I F ( K . LT . M ( I ) ) GO  TO  4 

370  o " =K— M ( I) 

380  4 K=;<+M(II) 

300  K=0 

400  DO  7 J=1  ,LX 

410  IF(K.LT.J)  GO  TO  5 

420  HOLD  = X(J) 

430  X ( J ) =X( K+ I ) 

440  X ( K+l  ) =HOLD 
450  5 DO  6 1=1  ,14 
400  1 1 =1 

470  IF(K.LT.M( I ))  GO  TO  7 

480  O K=K-M<  I) 

490  7 K =K+M  ( I I) 

500  IF ( SIGN. LT. 0.0)  RETURN 

510  DO  d 1=1 , LX 

520  8 X(  I)=X(  D/FLX 

530  RETURN 

540  END 


7H 


NOP MAG 


IOO  SUBROUTINE  NOR  MAG  ( LX,  X) 
I 10  DIMENSION  XCLX) 

I 2S  p=0.0 

I 30  DO  10  1=1  , LX 

140  I 0 8=AMAX  I ( A RS  ( X(  I))  ,R) 

14b  IE( P.E0.O.0) RETURN 

IbO  DO  20  1=1 , LX 

100  20  X ( I ) =X<  I)/B 

1 70  RETURN 

ISO  END 


7r 


OLD  PLO 


IOO  SUBROUT  INE  OLD  PLO(  U2 , F I , M ) 

I 10  DIMENSION  G2  ( M ) , F I (M) 

I I I DIMENSION  E( 04) ,0(64) 

119  CHARACTER  BLANK ( 63 )/63*"  "/ 

121  CHARACTER  STAR/"*"/ 

130  ROUND ( X) = X+ . 5 
132  DO  IS  I =1  ,M 
1 33  0( I)=G2(  I) 

134  IS  F ( I ) = E 1 ( I ) 

140  16  F0  = ( RIG< G, M) “SMALL ( G, M )) /42 . 

146  IF(FG.LE. 1 .0E-07)G0  TO  998 
ISO  DO  I K =1  ,M 
160  0(K)=G(K)/FG 
1 70  1 F(K) =2 . 0*F ( K ) 

180  CALL  GFSORTC  G, F ,M ) 

190  DO  19  K=1  ,5 
200  19  W R I i'E  ( 0 , 2 3 ) 

210  23  FORMAT ( 1 H ) 

220  01 =G< I ) *F0 
230  WRITE (0, I 00)01 
240  OLAST  =0 ( 1 ) 

2S0  DO  2 J=  l , M 

260  LODI F =ROUND( OLAST -G( J ) ) 

270  I F ( LGD  IF ) 8 , 8 , 6 

280  6 DO  7 L= 1 , LGD IF 

290  GG=( OLAST -FLOATC L) )*FG 

3 00  7 WRITE (0,  1 00)  GG 

310  100  FORMATUH  , E9.2,"  ") 

320  d LFD  IF=R()UND(  F(  J) +1  .) 

330  WRITE (0, 1 01 ) ( BLANKK  K) , K= 1 ,LFDI F) , STAR 
340  101  FORMATt 1H+,I0X,64AI ) 

35  t >. ' A 5T=0  ( J ) 

360  NR ITE ( 0, 2 00) 

370  2 00  FORM  AT(  I OX ,32 ( 7 I")) 

380  wRITE(  0,210) 

390  210  FORMATUOX,'  0123456789  11 

400  23  25  27  29  31  ') 

410  RETURN 

420  998  NR  ITE ( 0, 1 04 ) G(  1 ) 

430  104  FORMATC 'ALL  VALUES  ARE  EOUAL  T0',E9. 
440  RETURN 
500  END 


POP 


1 no  * I O 

1 in  FL 2 I N D , 0 
120  FL3IND.0 
I 30  *200 
140  ST A 
150  CLZE 
loO  CL. A 

I 70  i'AO  MT  ICKS 

180  CLAP 

190  CLA 

200  TAD  K5800 

210  CLOE 

220  CLA 

230  ST A 

240  DCA  FL2  TND 
250  STA 

260  DCA  FL3IND 
270  DCA  TALLY 
280  LOOP,  CLSK 
290  JMP  .-I 
295  CLSA 
300  JMS  SAMP 
310  ISZ  TALLY 
320  JMP  LOOP 
330  HLT 
340  SAMP,  0 
350  CLA 
360  TAD  K 2000 
3 70  JMS  ADCONV 
380  6221  / CDF2 
390  DCA  I FL2IND 
400  TAD  K2  00I 
410  JMS  ADCONV 
'°0  0231  /CDF 3 
4ou  DCA  I FL3IND 
440  Jf,.P  I SAMP 
450  ADCONV,  0 
460  A DSC 
470  A DCV 
480  ADSF 
490  JMP  .-1 
500  ADPP 

510  JMP  I ADCONV 
520  K200I  , 2001 
530  K2000,  2000 
540  TALLY , 0 
550  MTICKS,  -1415 
560  K5500 , 5500 
570  *300 
580  DISP,  LAS 


POP 

( contl 

590 

SMA 

0 00 

JMP 

. +3 

010 

6231 

/CDF3 

6 20 

JMP 

. +2 

630 

6221 

/CDF  2 

640 

CLA 

650 

DCA 

TALLY 

652 

TAD 

K2000 

654 

6552 

656 

CLA 

658 

6552 

660 

TAD 

I TALLY 

670 

0551 

680 

CLA 

690 

ISZ 

TALLY 

700 

JMP 

.-4 

710 

JMP 

DISP 

720 

$ 
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PLOT 


IOO  SUBROUTINE  PL0T(G2,F1 ,M) 

I 10  DIMENS  ION  G2(M  ) , FI  (M) 

111  DIMENSION  F( 64), 0(64) 

119  CHARACTER  BLANK ( 63 ) /63*"  «/ 

121  CHARACTER  STAR/"*"/ 

130  ROUND ( X)=X  + .5 
132  DO  15  1 = 1 , M 
1 33  G(  I) =G2( I ) 

134  15  F( I ) = F I ( I) 

140  16  FG=(BIG(G,M)-SMALL(G,M))/42. 

145  IF(FG.LE. I .0E-07)G0  TO  99B 
150  DO  I K =1  , M 
160  G(K)=G(K)/FG 
I 70  I F(  K) =2 . 0*F ( K ) 

180  CALL  GFS()RT(G,  F,M) 

I 00  DO  10  K = I ,5 
200  10  rtRITE( 0,23) 

210  23  FORMAT ( IH  ) 

220  Gl =G( I ) *FG 
230  rtRITECO, I 00)01 
240  GLAST=G(  I ) 

250  DO  2 J=l  ,M 

260  LOD I F=ROUND( GLAST-GC  J ) ) 

270  IF ( LODI F) 8, 8, 6 

280  6 DO  7 L=1  , LGD I F 

290  GG=( GLAST -FLOAT! L ) ) *FG 

300  7 rtRITECO, I 00)GG 

310  100  FORMAT ( 1 H , E9.2,"  ") 

320  8 LFDIF=R()UND(F(  J)  + l .) 

330  m R I TE ( 0 , 1 01  ) ( PLANK(K)  ,K=1  ,LFDI F)  f STAR 
340  101  FORMAT! 1 H+ , I 0X,64A 1 ) 

350  2 GLAST =0(J) 

360  rtRITE (0,200) 

■’■’O  2^0  FORMAT!  1 OX  , 32(  ' I')) 

380  4R  ITE(0, 21  0) 

390  2 1 0 FORMAT ( I OX,  ' 0 I 2 3 4 5 6 7 8 9 II  13  15  17  19  21'’, 
400  &'  23  25  27  29  31  ') 

410  RETURN 

420  998  «RITE( 0, 104)0( I ) 

4 30  104  FORMAT  ('ALL  VALUES  ARE  EQUAL  T0',E9.2) 

440  RETURN 
500  END 


PLOTTR 


1^0  SUBROUTINE  PLOTTR ( Yl , X,M) 

110  DIMENSION  Yl  (M)  ,X(  M)  ,Y  ( 1 28) 

120  CHARACTER  L INE ( 63)  /63*"  »/,RLAN/"  "/,  STAR/"*"/ , DOT/" . "/ 

129  DO  °o  1=1  ,M 

130  99  Y ( I ) =YI  ( I) 

131  SMALL=Y  ( 1 ) 

140  PIG=SMALL 

1 SO  DO  40  I =2  ,M 
loO  VALU  E= Y(  I ) 

170  IF(  VALUE-BIG)  20, 20, 10 
180  10  BIG= VALUE 
190  00  TO  40 

200  20  IF( VALUE-SMALL) 30,40,40 
210  30  SMALL= VALUE 
220  40  CONTINUE 

230  IF (ABS( PI G)-ABS( SMALL) >50,60,60 
240  50  SCALE =ABS( SMALD/31  . 

250  GO  TO  70 

260  60  SCALE =ARS( B IG)/3I  . 

270  70  I F ( BIG-SMALL)  I 00,80,100 
280  dO  hRITE! 0,90) BIG 

290  90  FOR:.iAT(  ' NO  GRAPH  CAN  BE  DRAWN  SINCE  ALL  VALUES  ARE  ',E15, 7) 

300  RETURN 

310  100  WRITE(0, I 10) PIG, SMALL, SCALE 

320  110  FORMAT  ('  LARGEST  VALUE  IS  ', El  5. 7,/, 'SMALLEST  VALUE  IS  ', 

330  &E15.7,/, / MULTIPLY  EACH  ORDINATE  READING  BY  THE  SCALE  FACTOR  ',EI5.7) 
340  DO  10  K=1  ,5 
350  10  HP  ITE(  0 , 23 ) 

360  23  FORMAT! I H ) 

3 70  DO  1 30  I = 1 , M 

380  130  Y ( I ) =Y!  I) /SCALE 
390  n,.*'E(0,2IO) 

400  210  FORMAT! 9X , '-30-27-24-2  1-18-15-12  -9  -6  -3  0 3 6 9 12', 

410  &'  15  18  21  24  27  30') 

420  NR  ITE(0, 1 20) 

430  120  FORMAT! 1 I X, 21! 'I  ')) 

440  DO  1 31  1=1 ,63 

450  I 31  LINE!  I ) =RLAN 
460  DO  1 0 00  1=1  , M 
470  VALUE  =Y(  I ) 

4 80  I ND EX  = 32 . + VALU  E 

481  INDEX=MAXO(  INDEX,  I ) 

482  INDEX  =MIN0(  INDEX, 63) 

490  LINE! 32)=D0T 

5 00  LINE!  I NDE  X)  =ST  AR 

510  WRITE (0, 1 80) X ( I) ,LINE 

520  LINE!  INDEX  )=PLAM 

530  180  FORMAT! F 7.2, 3X ,6 3A 1 ) 

540  1000  CONTINUE 
550  RETURN 


SO 


PLOTTR 


(continued) 


RB'.AV 

100  SUBROUTINE  REMAV(LY,Y) 
I 10  DIMENSION  Y(LY) 

1 27  S=U.O 
1 JO  DO  10  I =1  , LY 
1 40  1 0 S=S+Y( I ) 

ISO  A V=S/ELOAT( LY ) 
loO  DO  20  I =1 , LY 
170  20  Y ( I) =Y( I ) -A V 
100  RETURN 
190  END 


' 


F 


RZFRO 


100  SUBROUTINE  RZERO(LX.X) 
110  DIMENSION  X(  LX) 

I 20  DO  I 1 = 1 , LX 
130  I X(I)=0.0 
140  RETURN 
1 50  END 


SIMEQ 


14  00  SUBROUTINE  SI MEOC M, N , A , R, C) 

1410  * NMAX=LAROEST  VALUE  OF  N TO  RE  PROCESSED 
1420  * NONDUMMY  DIMENSION  S(NMAX,NMAX) 

1430  * FOR  EXAMPLE,  IF  NMAX=4  THEN 
1 440  DIMENSION  S(4,4) 

1 4b0  DIMENSION  A ( M , N ) , R( N , N ) , C ( M , N ) 

1 4oO  CALL  AO  VE(  N*N  , R, S ) 

1470  CALL  ESCAL( N,S, R) 

1480  DO  1 1=1  , M 

1 400  DO  1 J=1  ,N 
1 5 00  A ( I,J)  =0.0 
1510  DO  1 K=1  ,N 

1520  1 A( I ,J)=A(  I , J ) +C ( I ,K)*B(K,J) 

1530  CALL  ;AOVE(N*N,S,R) 

1540  RETURN 
1550  END 


I 


I 

S IMEOD 

loOO  SUBROUTINE  S IMEQU (M ,N , A, B, C,D ) 

1610  * NMAX=LARGEST  VALUE  OF  N TO  BE  PROCESSED 
1620  * NONDUMMY  DIMENSION  A(NMAX.NMAX) 

1630  * FOR  EXAMPLE,  IF  NMAX=4,  THEM 
1640  DIMENSION  S(4,4) 

1650  DIMENSION  A (M  , N ) , B(  N , II ) ,C(  M,  N ) 

1660  CALL  MOVE  ( N*N  , B , S ) 

1670  CALL  ESCALD(N,S, R,D) 

1 680  DO  1 1=1 ,M 
1 690  DO  1 J=1 ,N 
1 700  A ( I,J)=0.0 
1 710  DO  1 K=  1 , N 

1720  1 A( I ,J) =A(  I, J)+C( I ,J)*B(K, J) 

1730  CALL  MO  VE( N*N  ,S , B ) 

1 740  RETURN 
1 750  END 


SMALL 


I no  FUNCTION  SMALL! A , M) 
I 10  DIMENSION  A(M) 

1 20  S = A( I ) 

130  DO  I K =2 , M 
140  I F ( A ( K ) -S  ) 2 , 1 , I 
I SO  2 S=A ( K ) 

160  1 CONTINUE 
I o5  SMALL=S 
I 70  RETURN 
IBO  END 


' 


TRANS 


I 


100  SUBROUTINE  TRANS  ( P,  I D IMP,  Y , N , N V,  NR , L0G2NS ) 

I Ob  DIMENSION  Y(I024,NV) 

I 10  DIMENSION  P<  ID  IMP) 

; I 20  M=L0G2NS 

1 30  DO  I I =1 , NV-1 ,2 

140  I CALL  FASTI Y( 1 ,1 ) ,Y( J , I+l >,M,N) 

ISO  NT=N/2 

100  DO  8 1=1  tNV,2 

1 70  DO  8 J=2 , NT 
ISO  K=N-J+2 

190  Tl  =(  Y(  J,  I+l  ) -Y ( K , I + 1 )) 

200  T2=(Y(K,I )-Y(J,  I)) 

210  Y ( J , I ) =Y  ( J , I ) +Y ( K , I ) 

220  Y(J, I + l )=Y(J,I+1 )+Y(K,  I+l  ) 

230  Y ( K , I ) =T  1 
240  8 Y(K,  I+l  ) =T2 
2bO  DO  I 7 1 = 1 , N V 
260  Yd  , I ) =2 . 0*Y  ( 1 ,1  ) 

270  17  YCNT+1  , I ) =(  2.0)  *Y(  NT+I  ,I> 

280  NE=N/( 4*NR) 

2 90  1 1 M IN  =I4E+  1 
300  1 1 MAX=NT -NE 
310  I ISTE  P=NE+N  E 
320  L= I 

330  DO  2 I =1  , N V 
340  DO  2 J = I , N V 
350  P(L) =Y( 1 , I)  *YU  ,J) 

360  P( L+ 1 ) =0.0 
370  DO  3 K =2 , I IM  IN 
380  M=N -K+2 

390  3 P(L)=P(L>+2.0*( Y(K, I ) *Y ( K, J ) + Y( M, I)*Y(M,J)) 
400  L=L+2 

410  DO  4 II  = IIMIN,  UMAX,  1 1 STEP 
P(  L)  =0.0 
430  P(  L+ 1 ) 0.0 
440  KMAX=I  1+ IISTEP 
4bO  DO  5 K=I I , KMAX 
460  M=N-K+2 

4 70  P(L)  =P(L)+Y(K,  I )*Y(K,J)+Y(Mt  I)  *Y(  M,J) 

480  S PIL+1 )»P(L+l )+Y(K,I)*Y(M,J)-Y(M,I)*Y(K#J) 
490  4 l=L+2 

bOO  P(  L)  = Y(  NT+ 1 , 1 ) * Y(  NT+  I , J ) 

b I 0 P(L+I  ) =0.0 

b20  KM  IN=NT+ 1 -NE 

530  DO  6 K=KM  IN , NT 

540  M=N-K+2 

550  6 P(L)  = P(L)+2.0*(  Y(K,I  ) * Y(  K,  J) +Y  ( M , I ) * Y(  M,  J)  ) 
560  2 L=L+2 
570  RETURN 
580  END 


I 


UNPACK 


1 no  LET  1=0 

110  DIM  A(40O6) ,S( 13) 

120  FILE  #M  "PUB" 

130  FOR  J = 1 TO  LOEUI) 

140  READ  #1  : S$ 

150  CHANGE  S$  TO  S 
loO  FOR  K=  O TO  3 
170  LET  1=1+1 

180  LET  A.  ( I ) =1  6*S(  3*K+ 1 ) +1 NT( S( 3*K+2 ) / 1 6) 
190  LET  1=1+1 

200  LET  A ( I ) =256*M0U(S ( 3*K+2 ) , 1 6 )+S ( 3*K+3) 
210  NEXT  K 

2 20  NEXT  J 
230  SCRATCH  #1 
240  FOR  J=  1 TO  I 
250  m R ITE  #1  » A(  J) 

260  NEXT  J 

270  END 


PAR  Z 


100  SUBROUTINE  PARZ(M,W) 

110  DIMENSION  W ( M ) 

129  **  M IS  AN  EVEN  INTEGER 

1 30  DO  1 J=1 ,M/2  +1 

140  1 B ( J ) = I . -6 . *( FLOAT ( J-l ) /FLOAT CM)) **2+6 .*( FLOAT! J -I ) /FLOAT ( M ) )**3 
ISO  DO  2 J=M/2+2 ,M 

160  2 U ( J ) =2 .*( 1 . -FLOAT! J-l ) /FLOAT (M ) ) **3 
1 70  RETURN 
I 00  END 


ZERO 


I SUBROUTINE  ZERO ( LX, X) 
NO  DIMENSION  X(LX  ) 

I 40  COMPLEX  X 

130  I F( LX .LE.O) RETURN 

140  DO  I 1=1  , LX 

150  I X(I)=(0. 0,0.0) 

160  RETURN 
1 70  END 


LPSFLT1 


100  DIMENSION  rtT( 2005 ) , WTS( 49) , A I ( 2095 ) ,31 ( 2047 ) 
110  LIPkARY  "NI.OGN",  "CROSS",  "GENHTSM,  "FILTER" 

120  OPENFILE  I PUB",  "NUMERIC" 

I JO  R EAD( I ) < A I < J > ,J=I  ,2005) 

140  CALL  GENNTS(2,24,r<TS,49) 

IbO  CALL  F I LTER( 2095 ,49 , WTS, WT, A I , 204  7, B I , I 2 ) 

170  OPENF I LE  2,"NPUB", "NUMER IC" 

IBO  hR  ITE(  2)  Rl 
250  SLOP 
2oO  END 


LPSFLT2 


inn  DIMENSION  .ri  (2095)  ,WTS(49)  ,A2(  2005)  , R2(2047) 
110  LIBRARY  "NLOGN", "CROSS", "GENWTS", "FILTER" 

140  CALL  GENWTS ( 2, 24 , NTS ,40 ) 

190  OPENFILE  3, "CH 2 "."NUMERIC" 

200  READ(3)  < A2(  J ) , J =1  ,2095) 

210  CALL  FILTER<2095,49,KTS,WT,A2,2047,B2, 12) 

230  OPENFILE  4, "NCH2" , "NUMERIC" 

240  rt R IT E ( 4 ) B 2 
250  STOP 
260  END 


1 00 

'!.')!  )'  • ('(320  ) , !, 

1 (300)  ,0(  32)  , ; ( 33)  , r ( -1 

1 ! ) 

LI  * ■:  / "7  it " 

'7“ A V" , "ML 00 j “ , “CROSS" 

I -20 

>■  , »•  0‘  } L )T", 

11  2F50RT",  " I “SMALL *» 

1 30 

•'=32 

140 

I.Y=320 

1 b0 

OFF "FI  L : 3,  " FFOIJAT 

", "NUMERIC" 

IbO 

k-:ad(  3)  ( 

1 70 

CALL  ?:•'  ,2'V(LY,Y) 

1 JO 

CM  L C.  J M ( 370,  Y,  y 

,70,  if!  ,9) 

1 20 

CM  L.  ,<i>\RZ(  ",  7) 

700 

H-l ,/S4. 

>10 

• 1=  '+ 1 

>20 

"=(> 

?30 

cm  ' 7j';  :n  i',  • 

,21  , ),M|  ) 

> 10 

: ) 1 )- 1 , 23 

*>bO 

1 . ( 0)  = J-  1 

>40 

:\l  IJM1(  \r,  ') 

? JO 

■1  1 ( 0,  1 00)  (F(  J)  . 

0( J) , 3=1 , 33) 

7 / j 

1'X  !•  ) '.Arc  IH  >,F5. 

2,47, ! ; V. 2/) 

7J0 

00  7 J - 1 , / 0 

7 70 

2 . L A ;(  J)=j- 1 

100 

CALL  FLOTTRt  ! 1 , i’LA 

0,b0) 

110 

i.7i  l‘!:(  0 2 )0)(  i'L \ H 

J)  , f 1 ( J) , J=1 , /O) 

220 

?!)(’  1-0  ' i Vi'(  1 HO,.-'/. 

I ,47,  7-7.  >/) 

120 

pm  i , f 

240 

S 1'OP 

3 jO 

p in 

, )H(  33), TLA  H ■ 
..VARZ",  "St-FCi  " 
pLo  ~i 
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n LiU£-  .iOU  1 1 :!.£  RSFECT(H,*\N,R!  ) 

I 10  nr  u.!:3I  Oil  R 1 ( <1  ) ,0(  'I  ) 
i ?')  c(y,'LnX  c(ino) 
no  * "|=ih;  ?*o=?**n 

i i ; * ccm,  jmurc  ; jm  is  nowiJUmw  jimhmsio?:  >=  f.iicn 

IjQ  * '1  iiH  MAX  I MUM  I.  AO  M . 

I V)  C(  1 )=2Mi>LX(RI  ( I ) , 0 . ) 

I/O  DO  I J=2,M 

MO  I C(J)=L"U>LX(2.*HI(J),().) 

I VO  C(’M  )=CM3LX(I71  (Ml  ),0.) 

.200  no  2 J='+2,M+1 

on  2 c(  j)=co.  ,0. ) 

220  CALL  iJL  I .,-  + ■) 

030  no  3 j- 1 , 0 1 

240  3 r,(  J)=2.*H*FMSAL(C(  J ) ) 

250  RETtjf/ri 
260  E'lO 


6 bo  6 con  it; Ur: 

670  JO  7 1=1, UP  I 

/O')  L)()  / J =2 , 4 

vn  7 s(i,j,i)=s(i,i,j) 

167)  L)()  20  1 J = I , MV-  I 

160.)  DO  201  K=J  + 1 , NV 

16  10  DO  201  1 = 1,  Nil  1 

1620  CS5=CC.AB(S( I , J, J ) )*CCAB( S( I,K,K>> 

1630  IF( COS- I .OH-O / ) 17,10,10 
1640  17  S(r,K,J)=(0.0,0.0) 

1 6j0  00  ro  201 

1660  lo  S( I,K,J)=CMPLX(CCAB(5(I, J,K ) )**2/CSS,0.0) 

16/0  201  CONTINUE 

1600  DO  202  J 1 = 1 , NV-2 
1670  DO  202  J2=JI+ I , NV- I 
1700  DO  202  J 3= J2+ 1 , NV 

17  10  DO  202  1 = 1,  N!3  I 

1720  RES=REAL(S(  I,  J3,  J3  ) ) 

1730  SI  I 3 =S ( I , J 1 , J 1 )*( I .0-S ( I , J3, J 1 )) 

1740  S22  3=3  ( I , J2 , J2 ) * ( 1 .-S(  I , J3 , J2  ) ) 

1 760  I F ( RES- 1 . OE-O  7 ) 701,902,702 
1760  901  S I 23=S ( I, J1 ,J2) 

1770  GO  TO  703 

1780  702  S I 2 3=S ( I , Jl , J2)-S(  I , J I , J3 ) *CONJG(  3(  I , J2,  J3  ) )/RES 

1/70  703  I F ( CCAB( S 1 I 3 )*CCAB( S223 )- 1 .0E-07 ) 601,602,602 

1800  601  il?.(  I , J 1 , J2,  J3  ) =0.0 

1810  GO  TO  202 

1820  602  W2(  I , Jl , J2,J3)=(CCAB( SI  23 )**2 )/( CCAB( S 1 1 3 )*CCAB(S223 ) ) 

1830  RHS=REAL(S(I,J2,J2)) 

1840  SI I2=S( I, Jl , Jl )*( I .0-S( I , J2, Jl  )) 

I860  S332=S( I, J3,J3)*( l.-S( I,J3,J2)) 

I860  I F ( RES- I .OE-O 7 ) 1001,1002,1002 
13/0  1001  S I 32=S( I , J I , J3 ) 

1 330  GO  101003 

1870  1002  3I32=S(I,JI,J3)-S(I,JI,J2)*S(I,J2, J3)/RES 

1 700  1003  I F ( CCA3 ( S I I 2 ) *CCAB( S332 )- 1 . 0E-07 ) 701,702,702 

171)  701  ;V2  ( I , J I , J.3 , J2  ) =0 . 0 

1920  GO  i'O  202 

173  , 32 ( I , J 1 ,J3,J2)=( CCA3 ( S I 32 ) **2 )/( CCAB(S 1 I 2 )*CCAB( 5332 ) ) 

1740  PES=  EAL(S( I , Jl , Jl ) ) 

I960  S221=3( I, J2, J2>*( |.-S( I,J2,JI )) 

I960  S 33 1 =3 ( I, J3,J3)*( 1 .-S( I,J3,JI )) 

17/.)  Ir(  .3E5-  I . OE-O  / ) 1101,1102,1102 
I 980  1101  S23 1 =S( I , J2, J3 ) 

I 773  GO  TO  I 1 03 

2000  1 102  523  1 =S(  I,  J2 , J3  )-S  ( I , J I , J3  )*C()NJG(  S(  I , J I , J2  ) ) /RES 

2010  1103  I F ( CCAB (3221  )*CCA3 ( S 33 1 )- I . OE-O 7 ) 801,302,602 

2020  801  W2( I , J2, J3, Jl )=0.0 

20.30  GO  TO  202 

2040  802  )V2(  I,J2,J3,  Jl  ) = (CCAH(S23I  ) **2 ) /(  CCAB  C S22 1 )*CCA3(S.33I  )) 

2060  202  CONTINUE 

2060  CONTINUE 

2070  DO  71  J= 1 , NB I 

2080  71  F ( J)=J- I 

2090  DO  600  J= I , NV- 1 

2100  DO  500  K= J+ I , NV 

2110  A’R  I l‘E( 0, 300)CH(  J)  ,CH(  K ) 

2120  300  format ( * coherence  for  channels  ',ai,*-  and  ',ad 


98 


2130  DO  31 V 1 = 1, MB  I 
2 I 1)  SP(I)=REAL(S(I,K,J)) 

2 1-jJ  31-  CONTINUE 

3160  CALL  OLi)PT.C)(  SP,r , NB ) 

21/0  NR  I i 1:  ( 0,  I 03  ) 

2ldO  103  FORM AT ( 5 ( IH  ,/)) 

3100  DO  jOO  L=I,NV 
3300  I F ( < J-L ) *(  K-L  ) ) 4 | d , 500 , 4 1 d 
3310  4 Id  DO  b 12  1=1  , N'B  I 

3330  3P(  I ) = 12  ( I,  J,K,L) 

22.3  ) b I 2 CO’JTI  DUE 

2240  V;riTE(0,.30I  ) CM(  J ) , CH(  K ) ,CH(  L ) 

22.1*0  301  FORMAT  < * PARTIAL  COHERENCE  dETWEcN  CHANNELS  ',31,'  AND 
2260  *Al,/,'  AFTER  THE  INFLUENCE  OF  CHANNEL  ',AI,'  HAS  13  EE!  REMOVED') 
227  0 CALL  OLDPIJK  3P,  F,  MB  ) 

2230  r,RITE(0,  102) 

22>0  bOO  CONTINUE 
2 300  DO  417  J= I , NV 
2310  WEI TE(0, 3V7)CH( J) 

2320  317  FORMAT! ' AUTOSPECTRA  FOR  CHANNEL  ',AI) 

2330  I X=2*NH I *( NV I *( J— I )-((J-l )*(J-2)  )/?)-! 

2310  DO  4 16  1=  I tNI3 1 

23’j  > 4 16  JR(  I ) = Rl;AL(  S(  I , J,  J ) ) 

2 3oO  CALL  OLDPLO! SH ( I ) , F , N9 ) 

23  /0  4 17  JR  I Ti;(0,  1 0?  ) 

23rtD  STOP 
2 3 CO  HMD 


WMSSt M 


10)  dipemuio"  xu 000,2 ).z(  iooo) 

H ) l.IbrtARY  "rtNOI  SE" , "FOPTL!  b*** iFSAVFL" 
120  DATA  LZ/ I 000/ 

MO  00  I J=  1 , 2 

IV)  CALL  N *l() ISE(LZ,Z,B.0) 

160  DO  I 1=1,1000 
160  I X(I,J)=Z(I) 

I 6*j  CALL  FSAVrL(  "WNSDAT",  " ",  I ) 

I/O  OPE NF ILL  2, "WMSOAT", "NUMERIC" 

>60  NP  HE  ( 2 ) (<X(I,J),I  = I,I  000  ) , J=  I ,?) 
100  STOP 
200  FNL) 


/ 


/<70 


II  ) * tW  = N!J4H;R  Ur  CHANNELS 

I jo  * r r.-=,j!J"  HR  Or  FREOUFMCY  tiANDS(A  POWER  Oh  2) 

MO  * J iC  A ' * 5 IS  AL  JO  A POWER  OF  2 
IV)  * S . = SA  'FLI''0  R\TE=I/H 
160  Jr  X=INpJT  JEW  I E3(  ARRAY) 

I/)  * r=4RRA  { r Oft  S TORINO  CROSS  SPECTRA 
I n 01  •"  FUSION  <(  1024,2)  ,P(  I 9rt) ,S(  I 12)  ,C(  1.32  ) ,F(  33) 
lrj-j  UP  EMSl  00  5 1 ( 31,2,2) 

1 /)  CHARACTER  CH( 2) /" I ", "2"/ 

200  EON  I VALENCE  (5,0 

2)1  HOli  I VALENCE  (3,  51) 

21  ) DA'IA  MS.  MV,  NH.JSCANS,  SR,  PI/1000.  2,  32,  1024,64.,  1.  14159265/ 

2 IS  LI  PRAM Y "FAST", " TRAMS" , "MOVE" , "NOR MAO", "COHERE", "PLOT", “OFSOPT" 
2 I S A, "BIO", "SMALL" 

220  oPfcUFILE  2, "WMSOAT", "NUMERIC" 

210  REA1)(  2 ) ( ( X ( J,  I ) , J=  I , NS ) , 1 = 1 , MV) 

240  * DLTRH.NO  THE  SERIES  BY  SUBTRACTING  FROM  EACH  SERIES  ITS 
2‘jO  * I.LA51  SQUARES  LINEAR  REGRESSION  LIME 
260  ENS=NS 

2/0  Ti3AR=0.5*(rNS  + I .0) 

260  1 SUMS  )=( FM3* ( FNS+ I ,0)*(FMS+FNS  + l .0 ) )/6.0 

290  no  / 6 I2=|, MV 

100  SU**=0.0 

310  CRS>JRO=0 .0 

320  1)0  //  I 1 = 1 .MS 

310  S’JM  = SJM  + )((  11,12) 

140  77  CRSPRO  =CRSPRD+FLOAT(  I I )*X(  11,12) 
luO  FT-  AM=S.)M/FMS 

iso  FU-'i  A=(  CR S P R 0- r NS*T)3 A R*F *■< E A M ) / ( TSUMSQ-FMS*TI3A  R*TBAR) 

1/0  DO  70  I I =| , NS 

130  FRt:0=FM  EAN+PETA*(  FL()AT(  I I )-TB AR  ) 

190  7 3 X(  I I , 12 ) =X ( I I , 1 2)-FRE0 

400  76  C0M1IMUE 

410  * uiNOof  EACH  SERIES  WITH  A COSINE  TAPER 
420  IR=MS/10 
4 10  R = I R 

4 30  DO  / 9 ii  = i,il: 

4 )0  F I I = I 1 
4 "rj  FP,1=FI  1-0.5 

4 60  HI*  00  .'=0.3*(  I . O-COS ( PI  *t- 1 M T/R  ) ) 

3.  T3=NS+I-II 

4 30  |/1  oC  I 2=  I,  MV 

490  X ( I I , I2)=A.’IMi;0R*X(  11,12) 

500  bO  X ( 1 1 . I 2 ) =W  I Nl)OW*X  (11,12) 

•j  10  7 9 CONTINUE 
320  L002N3-) 

510  N SC  A NS  = I 

540  54  I F( MS.LE.NSCANS )00  TO  55 
550  L002NS=L002NS+I 
560  N S C A M S=  N SC  A N S + N S C A N S 
5 /O  00  i 0 54 

360  53  I F( NS. HO. NSC A MS )00  TO  74 

390  * ZERO  OUT  THE  LAST  PART  OF  EACH  SERIES  IF  THE  NUMBER  OF 

600  * SCAMS  IS  MOT  A POWER  OF  2 

610  I I P H0N=NS+ 1 

620  DO  75  I 1=1 IBEON.MSCAHS 

6 30  DO  75  I 2=1  ,NV 


10 1 


040  i‘j  x 1 1 1 , i2)=o.o 

630  74  CONiT  MJE 

650  I r ( ''()!)(  IV, 2))  /0,o?,  70 

670  * in  pn  or  SF!?  IES  IS  000,  FILL  A n,J,*”U  SERI  ES  WITH  ZEROS 

6:j  0 70  NVI=iV+l 

690  IK)  OS  i i = i,nsca:is 

703  OS  X ( I I , N V I )=0.0 

7 10  00  it)  j5 

720  OP  NVIaMV 

730  do  COM II NUE 

740  I [)I  MP=NV  I * ( N V I + 1 )*(Nb  + l ) 

7 50  CALL  fRANS ( P , I DI MP , X , MSCAMS,  MV  I , MO  , LOG  2 MS) 

/ 60  * CROSS  SPECTRUM  ESTIMATES  ARE  IN  ARRAY  P 

//O  * THE  CROSS  SPECTRUM  ESTIMATES  IN  ARRAY  P ARE  SCALED  BY 

/oO  * MULTIPLYING  BY  C 

7 VO  ANDPflR=EMS-l.25*R 

000  PSCANS=MSCAMS 
dIO  FNB=NB 

dPO  Er)=ESCAMS/(  FNI3+FN8 ) 

o30  C I =0 . 25/( SR* ( FD+ I ,0)*NNDPWR) 

040  I R0RSP  = ‘IB+NB+2 
350  ICOLSP=( NV  I *(  MV  I + I ) )/2 
060  I S IZEP  = I ROnISP*  I COLSP 
070  DO  95  11=1, ISIZEP 
OdO  95  P(  I I ) =C  I *P(  I I ) 

090  NO  I =MB+  I 
900  DO  1000  J= I , MVI 
910  00  1000  K=J , MV  1 
920  no  1000  1=1, NS  I 

930  SI  ( I , J,K)=P(2*NBI*(MV1*(  J- 1 )-('(  J- I )*(  J-2  ) J/2+K- J)  +1  +1-  l" ) 

9S-j  I F ( J-K ) 99,  I 000,  I 000 

940  99  SI  ( I J)=P(2*NBl*(!m*(  J-1>-(  ( .«■  * '*(  J-2  ) )/2+K-J)  +I  + I ) 

943  1000  CONTINUE 

950  CALL  COHEREC  NBI  ,VVI,S,C) 

960  no  / J = I , NB I 
970  7 E ( J ) = 3-  I 
9d0  Do  500  J=  I , NV-  1 

990  no  5 of)  :<=j+i,Mv 

1 00  ) MR  I Tfi  (0,300)  CH  ( J ) , CH ( ) 

ioi  > 300  format ( * coherence  for  channels  ',ai,'  and  lad 

10,.)  9,.:  PLOT(  C(  1 +NU I *( J- I ) + NB  I *NV  I * ( K-  I ) ) , F , NB  ) 

1043  100  rlRMATMH  ,F5 .0,  4X  , E9.  2 /) 

I 05  ) MR  ITEO.  102) 

1060  102  FORMA ;(5(!H  ,/)) 

10/  ) 500  CONTINUE 
lOd)  DO  105  J= I , NV 
lOd'j  UP  ITS (0,102) 

1093  WPITE(0, I07)CH( J)  " * 

1100  107  FORMA  T( ' AUTOSPECTRA  FOR  CHANNEL  ',A1) 

MI  3 CALL  PLOT ( C(  1 + NB  I *( J- 1 )+NR  I *NV  I *(  J-  I ) ) , F,  NB  ) 

1125  105  CONTINUE 
1 1 30  STOP 
1140  END 


IOZ 


WHOISTEST 


lf)0  * DIMENSION  X(N3*LX),RI  (LR*NS*NS),N!,M).F!M).S!  M*NS*NS) 

110  * DIMENSION  C(MI*NS*NS)  ,TLAG(2*LR-1  ) ,Z(  2*LR-  I ),  l'LAGA(LR) 

120  * NS=NUMBER  OF  CHANNELS 

130  * LX=NUMBER  OF  MTA  POINTS  FROM  EACH  CHANNEL 
140  * vi=M*l=LEMOTH  OF  TIME  LAC 
ISO  * LR=MAX I MUN  DESIRED  TIME  LAC  .LE.  LX 
160  * L=SMALLEST  INTEGER  SUCH  THAT  LX<=2**L 
I /O  * M=«/AXIMUM  LAG,  M=2**!  N- 1 ) 
no  * N IS  DEFINED  30  THAT  2*M=2**N 
I VO  * H “LENGTH  OF  SAMPLING  INTERVAL 
200  ★ L MX 3= LX* NS 
210  * LRMSNS=LR*NS*T> 

220  * V I NSNS=M I *N5*NS 
230  LIFE  NS  I ON  X(  1 000. 2 ) , R I ( 50,  2, 2 ) , H ( 32 ) ,F(  33), S(  142),  C(  33,2,2)  ,TLAG(  DD) 
240  *,Z< DD),TLAGA(50) 

250  CHARACTER  CH(  2)/'' I ",  "2"/ 

260  DATA  NS,  LX,  M,  Ml  , LR,L,  N,  LXNS,  LRNSNS , I N/2,  I 000 , 32 , 33 . 50 , 10 , 6, 2000 , 200,  2/ 
270  LIBRARY  "REMAV" , "NLOGN" , "CROSS", "WPARZ", "MACOR", "COQUAD", "MOVE" 

2B0  A, "NORM AC", "COHERE", "OLDPLO", "OFSORT", "BIG", "SMALL", "PLOTTR" 

2 DO  H= I ,/64. 

300  OPENFILE  2, "WNSDAT", "NUMERIC" 

310  READ C 2 ) X 

315  CALL  WPARZ! M , W ) 

320  DO  I J= I , NS 

330  1 CALL  REMAV( LX , X( I , J ) ) 

340  CALL  MACOR! NS, LX, X, LR, Rl , LXNS, LRNSNS, L) 

3-jO  CALL  COOUAD!  H,  NS,  M,  N,  W,  R I , S,  M I , LR) 

360  CALL  COHERE! Ml ,NS,S,C) 

370  DO  / J=  I , M I 

3 SO  7 F ( J)  = J-  I 

3 DO  DO  500  J= I , NS- I 

400  no  500  K=J+ 1 , NS 

A 10  WRITE(0,300)CH(  J),CH!K) 

420  300  FORMAT!  / COHERENCE  FOR  CHANNELS  ',A1,'  AND  LAI) 

4 30  CALL  OLDPLO!  C(  1,  J,K),F,M) 

450  100  FORMAT! IH  , 2 ( F6 . 2, 4X , ED. 2 , 6X ) /) 

460  rfRlTE(0, 102) 

170  102  FORMAT! 5 (1H  ,/)) 

4 BO  500  CONTINUE 
'♦  n DO  105  J=  I , NS 
500  WMTEtO,  I07)CH(  J) 

510  10/  FORMAT!'  AUTOSPECTRA  FOR  CHANNEL  ' , A I ) 

520  CALL  OLDPLO! C( 1, J, J),F,M) 

540  105  rtRI TEC 0,102) 

550  DO  700  1=1, NS- I 
560  DO  700  J= I + I , NS 
570  WRI I E(0, DO  I )CH( I ) , CH! J) 

-jdO  DO  I FORMAT!'  CROSS  CORRELATION  BETWEEN  CHANNELS  ',Al,'  AND  ' , A I ) 

5d0  DO  I OB  L=  I , LR-  1 

600  TLAG!  L)  =L-LR 

610  108  Z( L) =R I ( LR-L+I , J, I ) 

620  DO  IOD  L=LR,LR+LR- I 

6 30  TLAG!  L)=L-LR 

640  IOD  Z(L)*RI (L-LR+I , I , J) 

650  CALL  FLO ITR(Z ,TLAG, LR+LR- I ) 

660  WRITE! 0,102) 

6 10  700  CONTINUE 
660  DO  701  K= I , LR 


105 


6ro  7o i tlv;a<K)=k-i 

700  DO  bOI  J= I , NS 
/ I ) -NR  I i E(  0,  SO  I )CH(  J) 

/?.')  bO  I FORM  AT  ( / AUTO  CORRELATION  FOR  CHANNEL  ',AI) 
730  CALL  RLOffRC  Rl  ( I , J,  J)  , TLAGA.LR) 

740  dOI  WRI fE(0, I OR) 

/bO  PRINT, X 
760  STOP 
//J  EMU 


90  LIHkARY  "FILTER",  "SPKFLT" 

100  LI  PR  ARY  "SIMLA  T" , "OENFLT " , "ORMr'LI  " , "WKIOI SE",  "F0RTL1 P****FSAVFL" 

101  4,  "NLOG’l",  "CROSS " 

MO  * LZ.0E.LX*2*L.  L=50 
120  * LX.0E.LY  + 2*M.*I,  NW»5G 

130  DIMENSION  Z(  1000), X( 1000, 4), XF(o) ,YF(a) ,Y(500) 

140  & , A >V ( 101  ) ,0(  V)  ,,U(  1000)  , XZ(  1000,5)  ,C(4) 

145  A,A(4,t>),D(4,u),W(4,5) 

150  DATA  *3,  LZ , LX  , I.Y , WWT , 1.M2N/0 .015525,  1000,900,500.  101 , 10/ 

1 50  DATA  C( 1 ),C(2),C(3),C(4)/0.4,0. 1 ,0. 1,0. 1/ 

Is*)  DAI  A XF/4.  ,5.  ,7.  ,o.,  12..  13.,  15.,  16./ 

200  DAI  A YF/O., 1 ., I . , 0 , , 0 . , 0 .5 ,0 . 5,0. / 

210  DATA  A ( 1,  1 ),D(  1 , I ) ,/i(  I , 1 )/2.,0.,  I . / 

220  DATA  A(  I ,2) ,D(  I ,2) ,R(  1 ,2)/4 . ,6.,  I ./ 

2 30  DAI  A A ( I ,3).D( 1, 3 ) , W ( I, 3)/. 3.,  10.,  I ./ 

240  DAI  A 4(  1 ,4),D(  1 ,4) ,W( I ,4)/4.5, 14. ,2.  / 

250  DATA  A(  1 ,5)  ,.D(  1 .5)  ,W(  1 ,5)/4.  ,29.  , I ./ 

260  DATA  A(  2,  1 ) ,D(  2,  I ) ,.Y(  2,  1 )/.  1 , 10. , I ./ 

2/0  DAI  A A (2, 2),  i)  (2, 2),  (2, 2)/. 7,  20.,!./ 

2 SO  DAiA  A( 3, 1 ),  )( 3,  I ),R(3,  I )/.d, 3. , 1 ./ 

290  DATA  A(  3,2)  ,D(3,2)  , 'V(  3, 2 ) /.  2, 75  . . 1 ./ 

300  DATA  4(4.  I ),U(4,  I ),WM,  I )/,25,  12.,  I ./ 

310  DATA  A(4,2),D(4,2) ,rt(4,2)/0. , 1 . , 1 ./ 

320  DO  I K= I , 5 

340  1 CALL  SIM!)AT(H,A(  I ,K),D(  I ,K)  ,,<(  I ,K>  ,LZ,Z.LX,XZ(  I ,K  ) . W,  AW,  WT,  LN2N) 
345  CALL  'VNO I Sc*  ( LZ,  Z,  C ( 1 ) ) 

350  DO  2 J= 1 , LX 
360  X ( J , I ) =Z(  J ) 

370  DO  2 K=l  ,*j 
3 SO  2 X(  J,  I ) =X ( J, 1 )+XZ( J.K) 

390  CALL  ORMF LT( L Y , Y , NWT , AW , LZ , X ( 1 , I ), H, XF, YP ,8, 9, 0, AT. LN2H ) 

400  DO  3 N=2 , 4 
4 10  DO  4 < = I .2 

4 20  4 CALL  S I MDA i ( H , A ( N , K ),  D ( M , K ),  '.'(  N , K ),  LZ , Z , L Y , XZ ( I , K ) , NRT,  Art.WT,  LN?*J) 
430  CALL  AINOI SE(  LZ, Z, C(  M ) ) 

440  DO  3 J= I , LY 
450  X ( J, N ) =0 . 

46ij  0 5 K=l  ,2 

470  5 X( J,N)=X(J,N)+XZ( J,K) 

4 SO  3 X(J,M)=X(J,’J)+Z(J)  + Y(J) 

490  CALL  FSAVFLf  "FCHDAT",  " ",I) 

500  openfilf:  2,  "fchdat",  "numeric 

510  WRITE! 2) ( (X(I , J),I=1 ,LY) ,J=1 ,4) 

520  STOP 
530  END 


10  ) SUbrfOUi  I\'k  O.v Mr-L i ( LY.  Y.  .W,  AW.  LX,  X.H,  Xr,  YH ,U.U2,0. Wf , LM2N) 
I I ">  !)I»'i:N:iI  O J Y(  LY) , X(  LX)  , A.'H  ‘/WT  ) , MT(  I.X  ) ,Xr(  L.V ) , YE  ( Lft  ) ,0(  Lrt? ) 

I 20  * LX.0H.LY+2*N.V 
130  * LhP.GH.LV+l 
140  M.N  = ( fVi  T-  1 ) /2 
I-jO  N,n=‘LY+l 

160  CALL  GOMEL  PH.  Art  PM  I >,WltXF,YHtT  rt,0,L.<2) 

I/O  DO  I J= I , rM 

ISO  I Art(  Mrt  I - J ) =A  >'l  ( M.'il+J) 

CALL  E I L THn (I.X,  MWT , Ah , I.  i , X , L Y , Y,  LN2N  ) 

200  HE'l'URN 
210  END 


iO(o 


SiJl  <0Jir'4  SPTrLl  (M,  A I I , A) 

Di'^Voi  ri  a c : : . . i ) 

P=6.?d3l  it>3 

C=A  1/(1^.  73V2;>*H*./) 

IF(u.EO. 0.0)00  TO  POD 
A ( I )=(  ?*h*,-J)**2 
C 1 ( i)—;l ) 

C2=P*H*  ) 

C3=P*4*(  ')+,< ) 

DO  I J = 2 , iJr'j  I 
K=J-  I 

A ( J ) =-  ( I ./.<**  2 ) * ( OOS(CI*K)-2.*COS(C2*K)+COS(C3*K) ) 

1 CONTINUE 
DO  2 J=  I ,N  I] 

2 A ( J ) =C*A( J ) 

RETURN 

200  CI=P*H*f< 

A ( 1 ) =( 0 . j ) *C I **2 
DO  3 J=2,N4| 

K=J-  I 

3 A ( J ) = ( - I . 0 ) * ( I . /K**2) * ( COS( C I *K )— 1 .0) 

DO  4 J=I  , AM  I 

4 A CJ ) =C*A( J ) 

R 1: 1 U R N 

EMD 
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i n suBKourrui  stidatch.m  ,u,r, lz,z, lx, x, inr, a-v <1. ln2»d 

II)  DIMENSION  Z(LZ),X(LX),  AH(NrtT) , UT( LZ) 

120  * LZ  .OE . LX+2*N//=LX+NIVT-  I 
130  NW  = ( N-U  - I )/2 
1 4 0 MiVI  =M.V+  I 

10)  CALL  iPKEL  1 ( H, A 1 , l)t  W,  N«'<  I , A*/(  NW  I ) ) 

160  DC)  I J = I , N.  / 

I/O  1 Art(M/ll-J)=A/l(  M.N1+J) 

160  CALL  VNOlSEC  LZ.Z,  I .0) 

1 vO  CALL  E I L I'E .<  ( LZ,  N/II , AW , MT , Z , LX  , X , LN2N  ) 

200  RE'i  URN 
210  HMD 
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PROGRAM  OUTPUTS 
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AUTOS PECTRA  FOR  CHANNEL  I 
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EEG  c 1?fi5S  COKKeLATION  FUNCTION  OBf^IMfP  ®y  USING- 
MULSFECT 


auto  correlation  for  channel  I 

LARGEST  VALUE  IS  .2638448E+07 
SMALLEST  VALUE  IS  I6740R5E+07 

MULTIPLV  EACH  ORDINATE  READING  BY  THE  SCALE  FACTOR  .85H122E+05 


30.00 


AUTO  CORRELATION  FOR  CHANNEL  2 

LARGEST  VALUE  IS  .2606558E+07 

SMALLEST  VALUE  IS  1 673505 E+07 

MULTIPLY  EACH  ORDINATE  READING  BY  THE  SCALE  FACTOR 


•80B5670E+05 


COHERENCE  FOR  CHANNELS  I AND  2 
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PARTIAL  COHERENCE  BETWEEN  CHANNELS  I AND  2 
AFTER  THE  INFLUENCE  OF  CHANNEL  4 HAS  8EEN  REMOVED 
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PARTIAL  COHERENCE  BETWEEN  CHANNELS  l AND  3 
AFTER  THE  INFLUENCE  OF  CHANNEL  2 HAS  BEEN  REMOVED 
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COHERENCE  FOR  CHANNELS  I AND  4 
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PARTIAL  COHERENCE  BETWEEN  CHANNELS  1 AND  4 
AFTER  THE  INFLUENCE  OF  CHANNEL  3 HAS  BEEN  REMOVED 
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COHERENCE  FOR  CHANNELS  2 ANP  3 
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PARTIAL  COHERENCE  BETWEEN  CHANNELS  2 AND  3 
AFTER  THE  INFLUENCE  OF  CHANNEL  4 HAS  BEEN  REMOVED 
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COHERENCE  FOR  CHANNELS  2 AND  4 
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